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ABSTRACT 


The project, entitled "Computational Control of Flexible Aerospace Systems", was 
granted by NASA Langley Research Center (NAG-1-1436) started from January of 1992 and 
continued for three years. The main objective of this project is to establish a distributed parameter 
modeling technique for structural analysis, parameter estimation, vibration suppression and control 
synthesis of large flexible aerospace structures. 

The major focus of the first-year (1992) research was the distributed parameter modeling 
of large flexible aerospace structures with complex configurations. As an example, the Low- 
power Atmospheric Compensation Experiment (LACE) Satellite Model had been used as the 
testbed and physical object. The main accomplishments can be summarized as follows. A 
physical research LACE Model with a 10 by 10 feet supporting frame had been built at NASA 
LaRC. The classical modal test had been conducted on this model so that the natural frequencies 
of the model were measured which were used to compare with the analytical results obtained from 
the distributed parameter model of the LACE Model. A comprehensive dynamic formulation of 
the distributed parameter model for tether-beam-rigid-body assembled structures by using transfer 
matrix method had been systematically derived. Although the mathematical model must be 
derived, depending on the feature of different satellite structures, the methodology developed in 
the first-year research is definitely suitable for all tether-beam-rigid-body assembled structures. 
The software PDEMOD were used to find the natural frequencies of the system. The analytical 
values were comparable to the experimental results. A technical report entitled "Distributed 
Parameter Formulation of LACE Satellite Model by Using Transfer Matrix Method" had 
been submitted to NASA LaRC m . A summary of the report had been presented at the 9th 
VPI&SU Symposium on Dynamics & Control of Large Space Structures p] . The last two-year 
research was conducted on two phases. The first phase was the completion of the current version 
of PDEMOD Code and its related documentation. The second phase was to analyze the dynamic 
properties of a two-dimensional ground-based manipulator facility at the NASA Marshall Space 
Flight Center (MSFC) under various configurations, and to develop a methodology for vibration 
suppression of the end-effector by using distributed parameter modeling. 

This report concentrates on the research outputs produced in the last two years. The main 
accomplishments can be summarized as follows. A new version of the PDEMOD Code had been 
completed based on several incomplete versions Mr. Taylor was working on just before he died. 
The verification of the code had been conducted by comparing the results with those examples for 
which the exact theoretical solutions can be obtained. A summary of the theoretical background 
of the package along with the verification examples has been reported in a technical paper 
submitted to the Joint Applied Mechanics & Material Conference, ASME, Los Angeles, June 28- 
30, 1995 P] . Correspondingly, a brief USER'S MANUAL had been compiled, which mainly 
includes three parts: (1) Input data preparation, (2) Explanation of the Subroutines, (3) 
Specification of control variables. 



Meanwhile, a theoretical investigation of the NASA MSFC two-dimensional ground-based 
manipulator facility by using distributed parameter modeling technique has been conducted. A 
new mathematical treatment for dynamic analysis and control of large flexible manipulator systems 
has been conceived, which may provide a embryonic form of a more sophisticated mathematic 
model for future modified versions of the PDEMOD Codes. This research has been reported in 

two technical papers [4 5] . 


EXECUTIVE SUMMARY 

The common approach used in the analysis of structural dynamics and the interaction with 
the control synthesis of large flexible aerospace structures is die finite element method. Although 
the finite element method has been widely accepted, the significant limitations still exist. The 
elements used in the finite element method are usually void of dynamics, such as massless axial 
springs. The consequence is that hundreds and thousands of elements are needed to represent 
large flexible structures in order to acquire analytical accuracy. Because of the computational cost 
and numerical inaccuracies involved in generating solutions of large number of equations, there 
is a practical limit to the accuracy of finite element dynamic models. The high order of t e 
structural model requires an "order reduction" process before a control system can be designed. 
Seemingly unimportant modes can be inadvertently eliminated which prove later to be significant 
to control system performance and stability pj. 

Distributed parameter modeling is being seen to offer a viable alternative to the fmite 
element approach for modeling large flexible space structures. Distributed parameter models have 
the advantages of improved accuracy, reduced number of modal parameters, the avoidance o 
modal order reduction, and especially, the ability to represent the structural and control system 
dynamics in the same system of equations. Continuum models have been made of several flexible 
space structures, which include the Spacecraft Control Laboratory (SCOLE)^, Solar Array g t 
Experiment^,, NASA Mini-Mast Truss [9] , the Space Station Freedom [10J , and recently, the Low- 
power Atmospheric Compensation Experiment (LACE) Satellite Model P] . A computer software 
package aiming at performing structural analysis and control system synthesis had been initialized 
and is primarily completed for its basic functions. 

The software package PDEMOD was initialized by Dr. Lawrence W. Taylor, Jr. at NASA 
Langley Research Center during the middle of the 1980's. The first release of his work on 
PDEMOD package was in 1987 (1U2 , 13] . The initial interest of the package PDEMOD was to 
model the structural dynamics of general spacecraft configurations by using the distributed 
parameter approach. A system of partial differential equations is formulated and connected at 
their boundaries. The equations of motion for any number of rigid bodies are written in the 
frequency domain and in terms of the coefficients of the sinusoidal and hyperbolic functions which 
comprise the mode shapes. Distributed parameter models can, therefore, be generated for any 
three dimensional configurations describable by partial differential equations joined at their 
boundaries. The manual labor of generating such models is therefore avoided. 



Because of Dr. Taylor's sudden passing away, it becomes a urgent task to summarize and 
sift his research achievements, and make it available to the other researchers. With the NASA s 
support, a new version of the PDEMOD Code has been completed during the past two years based 
on several incomplete versions of Dr. Taylor. A summary of the theoretical background of the 
package along with the verification examples has been reported in a technical paper p] . 
Summarily, a complex large aerospace structure is considered as an assembly of flexible beam 
elements and rigid bodies. Each flexible beam element is represented by four independent partial 
differential equations which exhibit lateral bending in two axes, axial deformation, and torsion. 
A system of partial differential equations is then formulated and connected at the elements 
boundaries based on the compatibility conditions. The equations of motion for any number of 
rigid bodies are written in the frequency domain and in terms of the mode shape parameter 
coefficients. The deflections, forces and moments for both ends of a single beam element can be 
described in terms of the spatial derivatives of the solutions of the corresponding PDE's, further 
expressed in terms of the same set of mode shape parameter coefficients. Distributed parameter 
models can therefore be generated for any three-dimensional configurations describable by PDE s 
joined at their boundaries. An accomplishment of this task may provide an opportunity to more 
researchers to apply distributed parameter modeling techniques to a variety of aerospace 

structures. 

The verification of the Code has been conducted by comparing the results with those 
examples for which the exact theoretical solutions can be obtained. Four verified examples were 
included in the package: Example 1 - Bending of a Cantilevered Beam; Example 2 - Bending of 
a Clamped-Clamped Beam; Example 3 - Bending of a Cantilevered Beam with a Tip-Mass M; 
Example 4 - Torsion of a Cantilevered Beam. Correspondingly, a brief USER'S MANUAL, had 
been compiled, which mainly includes three parts: (1) Input data preparation, (2) Explanation of 
the Subroutines, (3) Specification of control variables. 

By investigating the potential of the distributed parameter modeling technique, Dr. Taylor 
and the other researchers expected that the PDEMOD may be further developed in the following 
aspects: (1) structural dynamics, modal frequencies and mode shapes; (2) parameter estimation 
of modal characteristics; (3) structural damping; (4) control system dynamics; and (5) design 
optimization. But, only the first of the program had been completed and included in the current 
package. To extend the functions of the current package, a massive research is being conducted, 
which suggests to modify the mathematical model and global system generating procedure^, to 
develop methodology for control synthesis [4 !4 15] . Instead of using the coefficients of the solution 
functions, the transfer matrix may be used ’finally, which provides a much more convenient way 
to describe the state-vector transition from one point of the structure to the other. 

These tentative ideas have been included in the second-phase research. A theoretical 
investigation of the NASA MSFC two-dimensional ground-based manipulator facility by using 
distributed parameter modeling technique has been conducted. The MSFC facility is planned to 
conduct research in the berthing operation and, in general, research into the control of multibody 
configurations that are loosely coupled with flexible manipulator linkages [16] . A new mathematical 
treatment for dynamic analysis and control of large flexible manipulator systems has been 



conceived, which may provide a embryonic form of a more sophisticated mathematical model for 
future modified versions of the PDEMOD Codes. This research has been reported in two 
technical papers [4 5] . 


ENCLOSURES 

The enclosures of this report are listed as follows, which represent the research 
accomplishments of this project. 

1. The PDEMOD Code and the computed results for four verified examples; 

2. USER'S MANUAL; 

3. Reference Paper 3: "PDEMOD - A Computer Program for Distributed Parameter Estimation 
of Flexible Aerospace Structures, Part I: Theory and Verification"; 

4. Reference Paper 4: "A Method of Superposing Rigid-Body Kinematics and Flexible Deflection 
for End-Effector Vibration Suppression of a Large Flexible Manipulator System". 


SUGGESTION TO FUTURE RESEARCH 

As mentioned before, the PDEMOD may be further developed in the following aspects: 
(1) structural dynamics, modal frequencies and mode shapes; (2) parameter estimation of modal 
characteristics; (3) structural damping; (4) control system dynamics; and (5) design optimization. 
To do so, the following two significant modifications will be necessary. First, the mathematical 
model supporting the current version of the PDEMOD Code must be modified. In current 
mathematical model, the equations of motion for any number of rigid bodies were written in terms 
of the mode shape parameter coefficients of the corresponding beam elements where the rigid 
bodies were attached. In other words, the state variables were chosen as the coefficients of the 
sinusoidal and hyperbolic functions which comprise the solutions of the corresponding PDE's. 
This model is not suitable for control synthesis. Instead of using the coefficients of the solution 
functions, the transfer matrix may be used, which provides a much more straightforward way to 
describe the state-vector transition from one point of the structure to the other, directly using 
deflections, slopes, forces and moments at both ends of individual beam element as state variables. 
Since deflections are the controlled variables and forces and moments are the motivating variables, 
choosing them as state variables are definitely more suitable for control purpose. 

Second, the package must include system identification and parameter estimation as an 
important part of the whole procedure. In general, the model used in distributed parameter 
analysis is actually an equivalent model of the real structure described by a set of PDE's. To keep 
the equivalency, a reasonable criterion to judge the equivalency must be properly set up first. The 
criterion used in current version of the PDEMOD was to keep the equivalency in the sense that 
the static stiffness were approximately equal between the distributed parameter model and the real 
structure. The equivalent parameters of the distributed parameter model, such as, mass, stiffness, 
radius of gyration, etc., were then determined based on this criterion. Further, the characteristic 



equation is solved to obtain the dynamic properties of the structure. However, the stiffness of a 
complex structure is usually related to the frequencies. The static equivalency can not guarantee 
dynamic equivalency in general, and error may arise. To correct the errors, current PDEMOD 
package adjusted the equivalent parameters based on some experimentally measured frequencies 
without specified mathematical algorithm. This correction is largely dependent on the user's 
experience, different users obtain different results at times. To overcome this deficiency, it is 
necessary to provide an algorithm to precisely define the dynamic equivalency, such as maximum 
likelihood estimator. 
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INTRODUCTION 


The software package PDEMOD was initialized by the late Lawrence W. Taylor, Jr. at 
NASA Langley Research Center during the middle of the 1980’s. The first release of his work on 
PDEMOD package was in 1987. PDEMOD was initially developed to model the structural 
dynamics of general spacecraft configurations by using the distributed parameter approach, which 
consists of three-dimensional network of flexible beams and rigid bodies. The building blocks 
from which three-dimensional configurations can be constructed consist of (1) beams, which have 
bending in two directions, torsion, and elongation degrees of freedom, and (2) rigid bodies, which 
are connected by any network of beam elements. The full six degrees of freedom are allowed at 
either end of the beam. Rigid bodies can be attached to the beam at any angle or body location. 
The modified Bemoulli-Euler beam equation is used to represent the bending and the wave 
equations for torsion and elongation. 

A system of partial differential equations (PDEs) is formulated and connected at the 
elements’ boundaries based on the compatibility conditions. The equations of motion for any 
number of rigid bodies are written in the frequency domain and in terms of the coefficients of the 
sinusoidal and hyperbolic functions which comprise the mode shapes. The force and moment 
vectors for both ends of a single beam element can be described in terms of the spatial derivatives 
of the solutions of the corresponding PDE’s. Distributed parameter models can therefore be 
generated for any three-dimensional configurations describable by PDEs joined at their 
boundaries. The manual labor of generating such models is therefore avoided. 

Because of Dr. Taylor’s sudden demise, it becomes an urgent task to summarize his 
research achievements and make them available to the other researchers. This work is being 
continued and this review may provide an opportunity to more researchers to apply distributed 
parameter modeling techniques to a variety of aerospace structures. The verification of the code 
has been conducted by comparing the results with those examples for which the exact theoretical 
solutions can be obtained, for instance, a simple beam with various boundary conditions, etc. 

The theoretical derivation of the formulation and the verification of the code have been 
included in Shen, et al.pj. This USER’S MANUAL concentrates on the explanation of the 
package itself. The package PDEMOD mainly includes three parts; (1) Input data, which 
specifies structural configuration, mechanical properties of the consisting beams and rigid bodies, 
and the natural frequency range one would like to search, etc., (2) Main body of the package, 
which conducts the calculation specified by the formulation developed in Shen, et al.pj. and (3) 
Subroutines necessary for completing the calculation. Since the formulation has been clearly 
described in Shen, et al.p], this Manual emphasizes the first and the third parts. It should please 
the user to know that it is necessary to read Shen, et al. t i], before reading this Manual. The user 
should begin by solving some of the example problems given in Shen, et alqij. The user should 
then proceed to work on their own complex problems. 

The USER’S MANUAL and Shen, et al.pj are the basic technical specification reference 
documents for the PDEMOD package. Although these two reference documents are sufficient 
for some users, other references, e.g. [3-14], will give more details of this package. It is our 
recommendation that users review as many of these references as possible to gain a thorough 
understanding of the PDEMOD package. 
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Neither the late Lawrence W. Taylor, Jr., nor his working colleagues, the compilers of this 
package and the authors of some of the related papers, assume any responsibility for the validity, 
accuracy, or applicability of any results. Users must verify their own results. 
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USER’S GUIDE 


A. SUBROUTINES: 


There are 17 subroutines in the package PDEMOD. The subroutines used m PDEMOD 
can be categorized as two types. The nine common-purpose subroutines (TYPE-1 
SUBROUTINES) are contained in a SUBROUTINE LIBRARY called “LODLIB.F , which is 
the selected part of the SUBROUTINE LIBRARY “SYSPAC” (SYStems analysis programs 
PACkage) at NASA Langley Research Center (LaRC)pj. The SYSPAC is a data base at LaRC 
for the purpose of making experimental data and a selection of analysis algorithms availab e o 
interested researchers studying aerodynamics, flight mechanics, structural dynamics and system 
identification. The other eight subroutines (TYPE-B SUBROUTINES) are not included in 
SYSPAC and are programmed specifically for PDEMOD. Their format is consistent with those 


oftheTYPE-I SUBROUTINES. , , a 

All vectors and matrices used in the subroutines are expressed in a vector form, and have a 
common format. The first four elements of each vector are respectively: (1) the number of rows, 
(2) the number of columns, (3) the total number of elements, and (4) the data time interval. s 
format allows matrix information to be readily accessible in programming data analysis 
procedures, for calling numerous subroutines, in printing and in plotting. The 17 subroutines in 
PDEMOD are described as follows. 


TYPE-I SUBROUTINES: 


1. SUBROUTINE ADD (P, A, Q, B, C) _ 

Description: Two compatible matrices (A and B) are multiplied by scalars (P 

respectively) and then added: C=P*A+Q*B. 


and Q, 


2. SUBROUTINE MULT (A, B, C) 

Description: Multiplies two matrices: C=A*B. 

3. SUBROUTINE MAKE (A, B) 

Description: Copies B in A: A=B. 

4. SUBROUTINE SET (A, II, JJ) 

Description: Creates a null matrix with II rows and JJ columns: A-[0]. 

5. SUBROUTINE SPIT (A, B) 

Description: Labels and lists a matrix. 

6. SUBROUTINE TRANS (A, B) t 

Description: Generates a matrix transpose: B=A . 

7. SUBROUTINE TILDA (A, B) 

Description: Forms the matrix equivalent of a cross product from a vector {A) 3 *i, 
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B = A = 


0 -a 3 a 2 

a 3 0 -a, 

-~ a 2 a \ 0 - 


8. SUBROUTINE JUXTV (A, B, C) . . 

Description- Combines by juxtaposition two compatible matrices in a vertical direction: 

fAl 


9. SUBROUTINE IDENT(A,II) 

Description: Forms an identity matrix A with dimension II. 


TYPE-n SUBROUTINES: 

10. SUBROUTINE AFORM (W, A, BODREF, CONFIG, L, P, Q, R, T, INRT, MASS, 
DUM, DUN, DUO, DUP) 

Description: Forms the system matrix A. (See Shen, et al.[ij, Eq.34) 

11. SUBROUTINE BODFORM (CONFIG, NBEAM, BODREF, NBODY) 

Description: Determines rigid bodies’ connectivity and the reference coordinate systems. 

12. SUBROUTINE PFORM (W, L, Z, E1X, EIY, EAZ, EIP, MPL, IPL, FZ, AGKX, 
AGKY, PF, PM) 

Description: Forms the matrices P F and P M (see Shen, et al.pj, Eqs.16 and 17). 

13. SUBROUTINE QFORM (W, L, Z, EIX, EIY, EAZ, EIP, MPL, IPL, FZ, AGKX, 
AGKY, QU, QS) 

Description: Forms the matrices Q„ and Q, (see Shen, et al.[ij, Eqs.14 and 15). 

14. SUBROUTINE PQFORM (W, L, E, P, Q, DUM, DUN) 

Description: Combines the matrices P and Q, respectively, by juxtaposition for multi-body, 

multi-beam system. 

15. SUBROUTINE UPPER (A, DETA) 

Description: Determines the determinant of the matrix A by transforming it as a upper 
triangular matrix. 

16. SUBROUTINE DIAG (A, SHAPE) 

Description: Diagrams the mode shapes. 

17. SUBROUTINE WSEARCH (W, DW, A, DETNEW, DETOLD, WI) 

Description: Search for the values of circular natural frequencies. 
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B. VARIABLES AND ARRAYS: 


There are a number of variables and arrays are used in PDEMOD. They are defined 

below: 

VARIABLES: 

NBEAM: Number of beams. 

NBODY: Number of rigid bodies. . . . 

NA=12*NBEAM, where the number 12 indicates that there are 12 unknown coefficients for each 

beam element (see Shen, et al.[ij, Eq.13). 

Li: Length of the ith beam. 

ElXi: Bending rigidity (EI X ) of the ith beam. 

EIYi: Bending rigidity (EI y ) of the ith beam. 

EEPi: Torsional rigidity (GI P ) of the ith beam. 

EAZi: Axial rigidity (EA Z ) of the ith beam. 

MPLi: Mass per length of the ith beam. 

IPLi: Mass moment of inertia per length of the ith beam. 

FZi: Initial tensile force for the ith tether element. 

AGKXi: Radius of gyration about x-axis of the cross-section area of the ith beam. 

AGKYi: Radius of gyration about y-axis of the cross-section area of the ith beam. 

MASSi: Mass of the ith rigid body. 

W: Circular natural frequency. 

WINC: Increment of the value of W for iteration. 

NIW: Specified number of iteration. 


ARRAYS: 

CONFIG (NBEAM, 3): Denotes the structural configuration: Column No.l=Beam Identification 
Number (ID); Column No.2=Inboard Body Identification Number; Column No.3-Outboard Body 
Identification Number. A negative sign prior to the numbers in columns 2 and 3 indicates 
which beam is used to define body axes. For example, a two-beam, three-body system is shown 
in Figure 1. The definitions of “Inboard Body” and “Outboard Body” are depicted at the nght- 
hand sides of the body ID Numbers. The 2 by 3 matrix CONFIG for this example is numbered as 

shown in the Tablel. 


Table 1 Example of the Matrix CONFIG 


Beam’s ID Inboard Body’s ID Outboard Body’s ID 


1 

-1 1 

-2 

2 

2 ^ 

-3 





Body 1 (Earth) — Inboard Body of Beam 1 
Figure 1 Two-Beam, Three-Body System 

The negative signs prior to the numbers located at (1, 2) and (1, 3) indicate that both bodies 1 and 
2 are defined in the beam l’s coordinate system, while the negative sign prior to the number 
located at (2, 3) indicates that body 3 is defined in the beam 2’s coordinate system. 

BODREF (NBODY, 2): Defines body’s connectivity and reference coordinate. The number of 
rows of BODREF equals to the number of bodies. Each row provides the connectivity 
information of the corresponding body, consecutively. Column No.l indicates the beam’s ID, the 
connected body uses this beam’s coordinate system as the reference system. Column No. 2 
indicates the mutual location between the body and the reference beam. If the body is the inboard 
body of that beam, the number equals to zero and if the body is the outboard body of that beam, 
the number equals to one. Note that BODREF is not a input array. All the numbers will be 
produced by calling SUBROUTINE BODFORM. In fact, CONFIG has provided all necessary 
information. 

Ril (3, 3): The eccentric matrix at the attachment point between the ith beam and its inboard 
body. 

RiO (3, 3): The eccentric matrix at the attachment point between the ith beam and its outboard 
body. 

R (4+18*NBEAM): The eccentric matrix containing all of matrices by juxtaposition: Ril and 
RiO, i=l, NBEAM. 

Ti (3, 3): The coordinate-transformation matrix. 

INRTi (3, 3): The mass-moment-of-inertia matrix of the ith body. 

INRT (4+9*NBODY): The mass-moment-of-inertia matrix containing all of INRTi matrices, 
i=l, NBODY, by juxtaposition. 

E (4+9*NBEAM): Beams’ mechanical-property matrix. The first four elements are used for the 
common purpose as mentioned before. Each consecutive nine numbers represent one beam’s 
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mechanical properties, respectively, that is: E(5)=EIX1, E(6)-EIY1, E(7)-EAZ1, E(8) EIP1, 
E(9)=MPL1, E(10)=IPL1, E(1 1)=FZ1, E(12)=AGKX1, E(13)=AGKY1, and so on. 


MASS (4+NBODY): Body’s mass matrix. 


p [4+(3*i2)*NBEAM]: Forms matrix P which consists of all the matrices Pf and Pm for ea ch 
beam element consecutively by juxtaposition (see Shen, et alp], Eqs.16 and 17). Expressed in 
matrix form, the matrix P is constructed as 


[P] = 


[Pf, ^3xi2 

[P Ml ] 3xl2 For Beam 1 

[P Fj ^3x12 

[PmJjx .2 For Beam 2 


Q [4+(3*12)*NBEAM]: Forms matrix Q which consists of all the matrices Q u and Q s for each 
beam element consecutively by juxtaposition (see Shen, et al.pj, Eqs.14 and 15). Expressed in 
matrix form, the matrix Q is constructed as 


[Q] = 


[Qu, ^3x12 

[Q h L .2 For 

[Q|l 2 ^3x12 
[QJ 3 x ,2 For 


Be ami 
Beam 2 


A [4+(12*NBEAM)*(12*NBEAM)]: Forms system matrix A which consists of all the matrices 
A f and A M for each beam element consecutively by juxtaposition, while the matrices A F and A M 
are constructed by the element matrices Pf, Pm, Qu, and Q* in the way indicated in Shen, et al.pj, 
Eq.33. 

DETOLD, and DETNEW [4+(12*NBEAM)*(12*NBEAM)]: Determine the determinant of 
the system matrix A by iteration. DETOLD is the previous one, while DETNEW the up-dated 
one. When the value of the determinant is small enough to be considered as zero, the 
corresponding natural frequency is found. 


C. CONTROL VARIABLES: 

NPROB: The problem number the user chooses to solve. If NPROB=l, then PROBLEM #1 is 
solved, and so forth. The current version of PDEMOD has encluded four verification examples: 
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EXAMPLE 1 - Bending of a Cantilevered Beam; EXAMPLE 2 - Bending of a Clamped-Clamped 
Beam; EXAMPLE 3 - Bending of a Cantilevered Beam with a Tip-Mass M; EXAMPLE 4 - 
Torsion of a Cantilevered Beam. User can add his own problems into the package and assign 
corresponding problem numbers. 

EFREQ: Index for conducting natural frequency analysis. If IFREQ=1, the natural frequency of 
the system will be computed. 

ISHP: Index for conducting mode shape analysis. If ISHP=1, the mode shape functions will be 
computed. 


D. SPECIAL NUMBERS: 

The following special numbers are used in the package to define some boundary 
conditions, null mass, or infinite mass, etc., so that the tedious modal reconstruction can be 
avoided. 

1. For a null mass, or an imaginary mass, MASSi, at the free end of a cantilevered beam, the 
package defines MASSi=0.9*10 8 . Correspondingly, the diagonal elements of its mass-moment- 
of-inertia matrix INRTi should be defined as the same number (0.9* 10 -8 ). 

2. For infinite mass, such as, the mass of the foundation of a cantilevered beam, MASSi, the 
package defines MASSi=999999999.0 Correspondingly, the diagonal elements of its mass- 
moment-of-inertia matrix INRTi should be defined as the same number (999999999.0). 

3. For restrictions of one-direction deflection, e.g., bending about x-axis, the package defines that 
the rigidity to resist the deflection in this direction approaches infinity, i.e., EIXi=999999999.0. 

4. For the elements except tether element, the package defines the initial tensile force FZi=0.0. 
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PDFMOD Program 


Computer-Printout Codes 



123 

543 

707 

111 

c 

C 


c 

C 

C 

C 

C 


PROGRAM PDEMODl 

INTEGER BODREF, CONFIG, OBODY %jrivo „. 

REAL INRTl , INRT2 , INRT3 , INRT4 , MASS1 , MASS2 , MASS3 , MASS4 , 

1L1 , L2 , L3 , MPL1 , MPL2 , MPL3 , MASS , 

2IPL1 , IPL2 , IPL3 , L , INRT , INRTl , INTRO 
DIMENSION DUM( 1300) , DUN(1300), DUO (1300) ,DUP(1300) , 

2 INRTl (13 ) , INRT2 ( 13 ) , INRT3 ( 13 ) , INRT4 ( 13 ) , 

3R1I (13 ) , RIO (13) ,R2I(13) , R20 (13 ) , R3I ( 13 ) , R30 (13 ) , 

4E (40) , A (2004) ,DETOLD(40) , CONFIG (5, 3) , 

5DETNEW (40 ) , BWl (8) ,BW2 (8) ,BW3 (8) , 

7PF (40) , PM (40) , R ( 94 ) f 

9MASSU4) !iNRt 169) ',BW(34) , TBODY ( 13 ) , TBEAM(13 ) - QU(40) , QS ( 40) 
DIMENSION RI ( 13 ) , RO (13 ) , INRTl (13 ) , INRTO (13 ) , P (436) , Q (436) , 

1TIT (13 ) , TRT ( 13 ) , BODREF ( 10 , 2 ) ,R1(13) ,CGI(7) ,CGTOT(7) , 

2TI(13) , SHAPE (2004) , XJI (2004) ,Dl(454) , 

3DELF ( 40 ) , DELC { 40 ) , 

^ FORMAT^ 2 8 H THE SYSTEM MATRIX WILL HAVE, 13 , 16HROWS AND COLUMNS) 
FORMAT (415, E15 . 5) 

FORMAT (7E14 .5) 

FORMAT (7110) 


USER CAN CONTROL WHICH PROGRAM FUNCTIONS ARE ACTIVE (=1), 
IFREQ=1 
ISHP=0 

USER MUST SELECT THE PROBLEM NUMBER DESIRED, i . e . NPROB=l , etc . 
PROBLEM #1 IS THE CANTILEVER BEAM BENDING; 

PROBLEM #2 IS THE CLAMPED-CLAMPED BEAM BENDING; 

PROBLEM #3 IS THE CANTILEVERED BEAM BENDING WITH A TIP_MASS; 
PROBLEM #4 IS THE CLAMPED-CLAMPED BEAM TORSION. 


AND NOT ( =0 ) 


NPROB=2 

if (nprob. eq. 1) goto 1001 
if (nprob. eq. 2 ) goto 1002 
if (nprob. eq. 3 ) goto 1003 
if (nprob. eq. 4) goto 1004 
c 

C-PROBLEM 1: ********** 


c 

1001 

C 

C 

C 

C 

c 

c 

c 

C 


140 
1 ^3 


C 


CONTINUE , , , . 

PROBLEM #1 is the CANTILEVERED BEAM BENDING (nbeam-1) . 

THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 
COL#l=BEAM I.D., COL#2-INBOARD BODY I.D., COL#3=OUTBOARD BODY 
NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES. 

For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 
Body2 (a Null-mass). Bodyl uses the Inboard end of Beaml to 
define its axes. Body2 uses the outboard end of Beaml to define 


its axes . 

CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3 )=-2 
NBEAM=1 

DO 140 IBEAM = 1 , NBEAM 

PRINT 777, CONFIG (IBEAM, 1) , CONFIG (IBEAM, 2) , CONFIG (IBEAM, 3) 
NA=12 *NBEAM 

CALL BODFORM ( CONFIG , NBEAM , BODREF , NBODY ) 

DO 163 I = 1, NBODY 

PRINT 777,1, BODREF (1,1), BODREF (1,2) 

CALL SET(WI, 100, 1) 

WI (1) =0 . 

INPUT FOR BEAM#1 , BODY#l INBOARD AND BODY#2 OUTBOARD 
Ll=130. 

EIX1=40000000 . 

EIY1=999999999 . 0 



EIP1=999999999 . 0 
MPL1=. 09556 
IPL1= .2907 
EAZ1=999999999 . 0 
FZ1=0 . 

- AGKX1=999999999 . 

AGKY1=999999999 . 

CALL SET (Rll , 3,1) 

Rll (5) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll ,DUM) 

CALL MAKE (Rll ,DUM) 

CALL SET (RIO, 3,1) 

RlO (5) =0.0 
RlO ( 6) =0 . 0 
RIO (7) =0.0 
CALL TILDA (RIO, DUM) 

CALL MAKE (RIO, DUM) 

CALL SET (Tl , 3 , 3 ) 

Tl (5) =1 . 

Tl (9) =1 . 

Tl (13 ) =1 . 

C INPUT FOR BODY # 1 , THE FIXED END, THE EARTH. 

MASS1=999999999 . 0 
CALL SET (INRTl ,3,3) 

INRT1 (5) =999999999.0 
INRTl (9) =999999999 . 0 
INRTl (13 ) =999999999.0 
PRINT 7 07, MASS 1 

C INPUT FOR BODY# 2 , THE NULL-MASS. 

MASS2=0 . 00000009 
CALL SET (INRT2 ,3,3) 

INRT2 (5) =0 . 00000009 
INRT2 (9) =0 . 00000009 
INRT2 (13) =0 . 00000009 
PRINT 7 07 , MASS2 
CALL SET ( E , NBEAM , 9 ) 

E ( 5 ) =EIXl 

E (6) =EIY1 

E (7 ) =EAZ1 

E ( 8 ) =EIP1 

E (9 ) =MPL1 

E (10) =IPL1 

E ( 11 ) =FZl 

E (12) =AGKX1 

E (13 ) =AGKY1 

CALL SET (L, NBEAM, 1) 

L (5) =L1 

CALL SET (MASS, NBODY, 1) 

MASS (5) =MASS1 
MASS (6) =MASS2 
CALL MAKE ( INRT , INRTl ) 

CALL JUXTV ( INRT , INRT 2 , INRT ) 

CALL MAKE (T, Tl ) 

CALL MAKE (R, Rll) 

CALL JUXTV (R, RIO, R) 

W=4 . 0 
WINC=1. 01 
NIW= 600 
GO TO 1000 

~cr 

c -PROBLEM 2: ********** 

c 

1002 CONTINUE 

C PROBLEM #2 is the CLAMPED_CLAMPED BEAM BENDING (nbeam=l) . 



THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 
COL#l=BEAM I.D., COL#2 -INBOARD BODY I.D., COL#3=OUTBOARD BODY 
NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES, 

For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 
Body2 (EARTH, EITHER) . Bodyl uses the Inboard end of Beaml to 
define its axes. Body2 uses the outboard end of Beaml to define 

its axes. 

CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3) =-2 
NBEAM=1 

DO 240 IBEAM = 1 , NBEAM , 

PRINT 777, CONFIG ( IBEAM, 1) , CONFIG (IBEAM, 2) , CONFIG ( IBEAM, 3) 

NA=12*NBEAM 

CALL BODFORM ( CONFIG , NBEAM , BODREF , NBODY ) 

DO 263 I = 1, NBODY 

PRINT 777,1, BODREF (1,1), BODREF (1,2) 

CALL SET (WI, 100,1) 

INPUT _ FOR BEAM#1, BODY#l INBOARD AND BODY# 2 OUTBOARD 
Ll=130 . 

EIX1=40000000 . 

EIY1=999 999999 . 0 
EIPl=999999999 . 0 
MPL1=. 09556 
IPL1= .2907 
EAZ1=999999999 . 0 
FZ1=0 . 

AGKX1=9999 . 

AGKY1=9 9 9999999 . 

CALL SET (Rll ,3,1) 

Rll (5) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll , DUM) 

CALL MAKE (Rll, DUM) 

CALL SET (RIO, 3,1) 

RIO (5) =0 . 0 
RIO (6) =0 . 0 
RIO (7 ) =0 . 0 
CALL TILDA (RIO, DUM) 

CALL MAKE (RIO, DUM) 

CALL SET (Tl , 3 , 3 ) 

Tl (5 ) =1 . 

Tl (9) =1 . 

Tl (13 ) =1 . 

INPUT FOR BODY # 1 , THE FIXED END, THE EARTH. 

MASS1=9 99999 999 . 0 
CALL SET (INRTl ,3,3) 

INRT1 (5) =9999999.0 
INRTl (9) =9999999.0 
INRTl (13) =9999999 .0 
PRINT 707 ,MASS1 

INPUT FOR BODY# 2 , THE NULL-MASS. 

MASS2=999999999 . 0 
CALL SET ( INRT2 ,3,3) 

INRT2 (5) =9999999 . 0 
INRT2 (9) =9 999999.0 
INRT2 (13 ) =9999999 . 0 
PRINT 707,MASS2 
CALL SET (E, NBEAM, 9) 

E ( 5 ) =EIXl 
E ( 6 ) =EIY1 
E (7) =EAZ1 
E (8) =EIPl 
E ( 9 ) =MPLl 



E ( 10 ) =IPLl 
E (11) =FZl 
E ( 12 ) =AGKXl 
E (13 ) =AGKY1 
CALL SET ( L , NBEAM , 1 ) 

L ( 5 ) =Ll 

CALL SET (MASS, NBODY, 1) 

MASS ( 5 ) =MASS1 
MASS ( 6 ) =MASS2 
CALL MAKE (INRT, INRTl) 

CALL JUXTV ( INRT , INRT2 , INRT ) 
CALL MAKE (T, Tl) 

CALL MAKE (R, Rll) 

CALL JUXTV (R, RIO, R) 

W=25 . 0 
WINC=1 . 01 
NIW=300 
GO TO 1000 


c-PROBLEM 3 ********** 


c 

1003 

C 

c 

C 

c 

c 

c 

c 

c 

C 


340 

363 

C 


CONTINUE , . , , . . . 

PROBLEM #3 is the CANTILEVERED BEAM BENDING (nbeam=l) with a 

TIP-BODY CONNECTED. 

THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 
COL#l=BEAM I.D., COL#2 -INBOARD BODY I.D., COL#3=OUTBOARD BODY 
NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES 

For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 
Body2 . Bodyl uses the Inboard end of Beaml to 

define its axes . Body2 uses the outboard end of Beaml to define 


its axes. 

CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3) =-2 
NBEAM=1 

DO 340 IBEAM = 1, NBEAM , , 

PRINT 777, CONFIG (IBEAM, 1) , CONFIG (IBEAM, 2) , CONFIG (IBEAM, 3 ) 

NA=12 *NBEAM 

CALL BODFORM (CONFIG, NBEAM, BODREF, NBODY) 

DO 363 I = 1, NBODY 

PRINT 777,1, BODREF (1,1), BODREF (1,2) 

CALL SET (WI, 100,1) 

WI (1) =0 . 

INPUT FOR BEAM#1, BODY#l INBOARD AND BODY#2 OUTBOARD 


Ll=3 . 077 


EIX1=175 . 9644 
EIY1=999999999 . 0 
EIP1=999999999 . 0 
MPL1=0 . 012 


IPL1=1 . OelO 


EAZ1=999999999 . 0 


FZ1=0 . 

AGKX1=999999999 . 
AGKY1= 9 99999999 . 
CALL SET (Rll, 3,1) 
Rll (5) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll , DUM) 
CALL MAKE (Rll, DUM) 
CALL SET (RIO, 3,1) 
RIO (5) =0 . 0 
RIO ( 6) =0 . 0 
RIO (7 ) =-0.07 
CALL TILDA (RIO, DUM) 
CALL MAKE (RIO, DUM) 



CALL SET (Tl ,3,3) 

Tl (5) =1 . 

Tl (9 ) =1 . 

Tl ( 13 ) =1 

INPUT FOR BODY#l , THE CLAMPED_END , THE EARTH. 
MASS1=99 9999 999 . 0 
CALL SET ( INRTl ,3,3) 

INRTl (5) =999999999.0 
INRTl (9) =999999999.0 
INRTl <13 ) =9 99999999.0 
PRINT 707 ,MASS1 

INPUT FOR BODY# 2 , THE REFLECTOR. 

MASS2=4. 952/32 .2 
CALL SET ( INRT2 ,3,3) 

INRT2 (5) =2 . 341e-4 
INRT2 (9 ) =2 . 341e-4 
INRT2 (13 ) =1 . 524e-3 
PRINT 707 ,MASS2 
CALL SET ( E , NBEAM , 9 ) 

E (5) =EIX1 

E ( 6 ) =EIY1 

E (7 ) =EAZl 

E ( 8 ) =EIPl 

E (9 ) =MPL1 

E ( 10 ) =IPLl 

E (11) =FZl 

E (12 ) =AGKXl 

E ( 13 ) =AGKY1 

CALL SET (L, NBEAM, 1) 

L (5) =Ll 

CALL SET (MASS , NBODY , 1 ) 

MASS (5) =MASS1 
MASS ( 6 ) =MASS2 
CALL MAKE (INRT, INRTl) 

CALL JUXTV ( INRT, INRT2 , INRT) 

CALL MAKE (T, Tl) 

CALL MAKE ( R , Rl I ) 

CALL JUXTV ( R , RlO , R ) 

W=10 . 0 
winc=l . 01 
niw=150 
GO TO 1000 


c-PROBLEM 4 ********** 


c 

1004 

C 

c 

C 

C 

C 

c 

c 

c 

C 


463 


PROBLEM #4 is the TORSION Of A CLAMOED-CLAMPED BEAM (nbeam 1) 
WITHOUT BODY CONNECTED. 

THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 
COL#l=BEAM I.D. COL#2 -INBOARD BODY I.D., COL#3— OUTBOARD BODY 
NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES. 

For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 
Body2 (another clamped end). Bodyl uses the Inboard end of Beaml to 
define its axes. Body2 uses the outboard end of Beaml to define 
its axes . 

CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3 )=-2 
NBEAM=1 

DO 440 IBEAM = 1, NBEAM 

PRINT 777 , CONFIG (IBEAM, 1) , CONFIG ( IBEAM, 2 ) , CONFIG ( IBEAM, 3) 
NA=12*NBEAM 

CALL BODFORM ( CONF IG , NBEAM , BODREF , NBODY ) 

do 463 i=l , nbody t „ 

print 777, i,bodref (l, 1) ,bodref (l, 2) 

CALL SET(WI, 100,1) 



WI (1) =0 . 

INPUT FOR BEAM 1, BODY#l INBOARD, BODY#2 OUTBOARD 
Ll=130 . 

EIX1=999999999 . 0 
EIyl=9 99999999 . 0 
EIpl=4 00000000 . 0 
MPL1=0 . 09556 
IPL1= .2907 
EAZ1=999999999 .0 
FZ1=0 . 

AGKX1=999999999 . 0 
AGKyl=999999999 . 0 
CALL SET (Rll/ 3,1) 

Rll (5) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll , DUM) 

CALL MAKE (Rll, DUM) 

CALL SET (RIO, 3,1) 

RIO (5) =0 . 

RIO (6) =0 . 

RIO (7 ) =0 . 

CALL TILDA (RIO, DUM) 

CALL MAKE (RIO, DUM) 

CALL SET (Tl ,3,3) 

Tl (5) =1 . 

Tl (9 ) =1 . 

Tl (13 ) =1 . 

INPUT FOR BODY l(one clamped end) 

MASS1=999999999 . 

CALL SET ( INRTl ,3,3) 

INRT1 (5) =999999999 . 0 
INRTl (9) =999999999 . 0 
INRTl (13 ) =999999999.0 
PRINT 707, MASS1 

INPUT FOR BODY 2 (another clamped end) 
MASS2=999999999 . 0 
CALL SET ( INRT2 ,3,3) 

INRT2 (5) =999999999.0 
INRT2 (9) =999999999.0 
INRT2 (13) =999999999.0 
PRINT 707, MASS2 
CALL SET ( E , NBEAM , 9 ) 

E ( 5 ) =EIXl 

E ( 6 ) =EIY1 

E (7 ) =EAZl 

E (8) =EIP1 

E (9) =MPL1 

E ( 10 ) =IPLl 

E (11) =FZl 

E (12) =AGKX1 

E ( 13 ) =AGKYl 

CALL SET (L, NBEAM, 1) 

L ( 5 ) =Ll 

CALL SET (MASS , NBODY , 1 ) 

MASS ( 5 ) =MASS1 
MASS (6) =MASS2 
CALL MAKE (INRT, INRTl) 

CALL JUXTV ( INRT , INRT2 , INRT ) 

CALL MAKE (T, Tl) 

CALL MAKE (R, Rll) 

CALL JUXTV (R, RIO, R) 

W=280 . 0 
WINC=1. 01 
NIW=300 
GO TO 1000 



***** TH e COMMON PART OF THE MAIN PROGRAM 


***** 


c 

c 

c 


1 ^00 


C 

C 


1 


31 


92 


30 

29 

28 

3 

999 

998 


CALL SPIT (L, 2H L) 

CALL SPIT (E, 2H E) 

CALL SPIT (R, 2H R) 

CALL SPIT (T, 2H T) 

CALL SPIT (MASS, 5H MASS) 
CALL SPIT ( INRT , 5H INRT) 
CALL SET ( DETNEW , NA , 1 ) 
DO 1 IW=1 , NIW 
DW=W* (WINC-1 . ) 


W=W+DW 

FORM FORCE AND MOMENT MATRICES "PF" AND "PM" 

FORM LINEAR AND ANGULAR DEFLECTION MATRICES "QU" AND 
CALL PQFORM (W, L, E, P, Q, DUM, DUN) 

CALL SET (A, NA, NA) 

CALL AFORM (W, A, BODREF , CONFIG , L, P, Q, R, T, INRT , MASS , 
1DUM , DUN , DUO , DUP ) 

CALL ADD{1 . , A, 0 . , A, A) 

CALL WSEARCH ( W , DW , A , DETNEW , DETOLD , WI ) 


CONTINUE 

AD=l./6. 283185 

CALL ADD ( AD , WI , 0 . , WI , F I ) 

CALL SPIT (FI, 3H FI) 

IF(ISHP)998, 998,31 
CONTINUE 
NW=WI (1) +. 01 
DO 3 IW=1 , NW 
W=WI (IW+4) 

CALL PQFORM (W, L, E, P, Q, DUM, DUN) 

CALL AFORM (W, A, BODREF , CONFIG , L, P, Q, R, T, INRT , MASS , 
1DUM , DUN , DUO , DUP ) 

CALL UPPER (A, DETNEW) 

DETN=1 . 

DO 92 I = 1,NA 
DETN=DETN* DETNEW ( 1+4 ) 

CALL SET ( DUM , 1 , NA ) 

DUM(NA+4) =1 . 

CALL DIAG ( A , DUM ) 

IF ( IW- 1)29,30,29 
CALL MAKE (SHAPE, DUM) 

GO TO 28 

CALL JUXTV( SHAPE, DUM, SHAPE) 

CONTINUE 

CONTINUE 

IF(ISHP) 998, 998, 999 
CALL SPIT (SHAPE, 5HSHAPE) 

STOP 

END 


" QS " 


c ***** ***** ***** 

c SUBROUTINES 

c ***** ***** ***** 

SUBROUTINE AFORM (W, A, BODREF, CONFIG, L, P, Q, R, T, INRT, MASS , 
1DUM, DUN, DUO, DUP) 

INTEGER APPEND , BEAMl , BODREF , CONFIG , OBODY 
REAL INRT, INRT I, INRTO, MASS, L 

DIMENSION DUM (13 00) , DUN (13 00) , DUO (13 00) , DUP (13 00) , 

2 INRT (69) , INRTI ( 13 ) , MASS (14) , L (9) , R (94) ,RI (13 ) , Rl (13 ) , 
3A(2004) , P (43 6 ) ,PF(40) ,PM(40) ,Q(436) , QU (40) , QS (40 ) , 

~ 4BODREF (10,2) ,CONFIG(5,3) ,T(49) , TBEAM ( 13 ) ,T1(13) , 

5TIT ( 13 ) , TRT ( 13 ) ,QUl(40) ,QSl(40) , TI (13 ) , RO (13 ) , 

6 INRTO (13 ) 

707 FORMAT (7 El 5 . 4 ) 

NBODY =MAS S ( 1 ) + . 0 1 



NBEAM=L ( 1 ) + . 0 1 
NA=NBEAM* 1 2 
APPEND=NBODY 
DO 2 IBEAM = 1 , NBEAM 
CALL SET (TBEAM, 3,3) 

— CALL SET (RI , 3 , 3 ) 

CALL SET ( RO ,3,3) 

DO 13 IBLOCK =1,9 

TBEAM ( IBLOCK+4 ) =T ( (IBEAM*3-3) *3+IBLOCK+4) 

RI ( IBLOCK+4 ) =R ( (IBEAM*6-6) *3+IBLOCK+4) 

13 RO( IBLOCK+4 )=R( ( IBEAM*6-3 ) *3+IBLOCK+4 ) 

C FORM MATRIX "A” FOR INBOARD BODY 

CALL SET(QU, 3,12) 

CALL SET(QS, 3,12) 

DO 12 IBLOCK = 1,36 

QU (IBLOCK+4 ) =Q ( ( IBEAM* 12- 12 ) *12+ IBLOCK+4 ) 

12 QS (IBLOCK+4) =Q ( (IBEAM*12-9) *12+IBLOCK+4 ) 

CALL TRANS (TBEAM, TIT) 

IBODY=CONFIG ( IBEAM, 2 ) 

IF(IBODY) 6, 4, 4 
6 IBODY=-IBODY 

ENREF=1 . 

GO TO 18 
4 ENREF=0 . 

BEAMl =BODREF ( I BODY, 1 ) 

IO=BODREF ( I BODY, 2 ) 

IF (BEAMl -IBEAM) 70,71,70 

70 CALL SET (Tl ,3,3) 

CALL SET (QUl ,3,12) 

CALL SET (QS1 ,3,12) 

CALL SET (Rl , 3 , 3 ) 

DO 72 IBLOCK =1,9 

Tl (IBLOCK+4) =T( (BEAMl*3-3) *3+IBLOCK+4) 

Rl (IBLOCK+4) =R( (BEAMl*6+3 *10-6 ) *3 + IBLOCK+4) 

DO 73 IBLOCK =1,36 

QUl (IBLOCK+4 ) =Q( (BEAMl*12+6*IO-12 ) *12+IBLOCK+4) 

73 QSl (IBLOCK+4) =Q( (BEAMl *12 +6 *10- 9 ) *12+ IBLOCK+4) 
CALL MULT (Tl, QUl, DUN) 

CALL ADD ( 1 . , RI , -1 . , Rl , DUM) 

CALL MULT ( DUM, Tl, DUO) 

CALL MULT (DUO, QSl, DUM) 

CALL ADD ( - 1 . , DUN , - 1 . , DUM , DUO ) 

CALL MULT (TIT, DUO, DUP) 

CALL MULT (Tl, QSl, DUM) 

CALL MULT (TIT, DUM, DUN) 

CALL ADD ( - 1 . , DUN , 0 . , DUN , DUN ) 

CALL JUXTV ( DUP , DUN , DUO ) 

CALL JUXTV ( QU , QS , DUP ) 

DO 75 IBLOCK =1,6 
DO 74 JBLOCK =1,12 

ID= (APPEND* 6+ IBLOCK- 1) *NA+JBLOCK+BEAMl*12-8 

74 A ( ID) =DUO ( (IBLOCK-1) *12+JBLOCK+4 ) 

DO 75 JBLOCK =1,12 

ID= (APPEND* 6+ IBLOCK-1) *NA+JBLOCK+IBEAM*12-8 

75 A ( ID) =DUP ( (IBLOCK-1) *12+JBLOCK+4 ) 
APPEND=APPEND+ 1 

71 CONTINUE 

18 CONTINUE 

EM=MASS ( IBODY+4 ) 

CALL SET ( INRTI ,3,3) 

DO 14 IBLOCK = 1,9 

IT4 INRTI (IBLOCK+4) =INRT( (IBODY*3-3) *3 + IBLOCK+4) 

CALL SET (PF, 3, 12) 

CALL SET (PM, 3, 12) 

DO 31 IBLOCK =1,36 

PF (IBLOCK+4) =P ( ( IBEAM* 12 -12 ) *12+IBLOCK+4 ) 



31 PM ( IBLOCK+4 ) =P ( (IBEAM*12-9) *12+IBLOCK+4) 

W2=W*W 
W2IN=1. /W2 

CALL MULT (RI , TBEAM, DUM) 

CALL MULT (TIT, DUM, TRT) 
tr' FIRST BLOCK, FORCE EQUATIONS, INBOARD 
IF (EM-999999999 .) 86, 87,86 

87 AD=ENREF 
AE=0 . 

GO TO 88 

86 AD=ENREF*EM 

AE=W2 IN 

88 CONTINUE 

CALL ADD (AD, QU, AE, PF, DUO) 

CALL MULT (TRT, QS, DUN) 

CALL ADD ( 1 . , DUO , +AD , DUN , DUO ) 

C SECOND BLOCK, MOMENT EQUATIONS, INBOARD 

CALL MULT (RI, TBEAM, DUM) 

CALL MULT (DUM, PF, DUN) 

CALL MULT (TBEAM, PM, DUM) 

CALL ADD ( +W2 IN , DUM , +W2 IN , DUN , DUP ) 

CALL MULT (TBEAM, QS, DUM) 

CALL MULT (INRTI, DUM, DUN) 

CALL ADD ( 1 . , DUP , ENREF , DUN , DUP ) 

CALL JUXTV ( DUO , DUP , DUO ) 

DO 5 IBLOCK =1,6 

DO 5 JBLOCK =1,12 . 

5 A ( ( IBODY*6+IBLOCK-7 ) *NA+IBEAM*12+ JBLOCK-8 ) =DUO ( ( IBLOCK- 1) 

1*12 + JBLOCK+ 4 ) 

C FORM MATRIX "A" FOR OUTBOARD BODY 

CALL SET (QU, 3,12) 

CALL SET (QS, 3,12) 

DO 19 IBLOCK = 1,36 

QU(IBLOCK+4) =Q( (IBEAM*12-6) *12 + IBLOCK+4) 

19 QS ( IBLOCK+4 ) =Q ( (IBEAM*12-3 ) *12+IBLOCK+4) 

OBODY=CONFIG ( IBEAM, 3 ) 

IF (OBODY) 15,15,16 

15 ENREF =1 . 

OBODY= -OBODY 
GO TO 17 

16 ENREF =0 . 

BEAMl =BODREF (OBODY, 1) 

IO=BODREF ( OBODY , 2 ) 

IF (BEAMl -IBEAM) 80,81,80 
80 CALL SET (Tl , 3 , 3 ) 

CALL SET (QUl ,3,12) 

CALL SET(QS1, 3, 12) 

CALL SET (Rl , 3 , 3 ) 

DO 82 IBLOCK = 1,9 

Tl (IBLOCK+4) =T ( (BEAMl*3-3) *3+IBLOCK+4) 

82 Rl( IBLOCK+4 )=R( ( BEAMl * 3 - 3 ) * 3 + IBLOCK+ 4 ) 

DO 83 IBLOCK =1,36 

QUl (IBLOCK+4) =Q( (BEAMl *12 +6 * 10-12 ) *12+IBLOCK+4) 

83 QS1 (IBLOCK+4) =Q( (BEAMl*12+6*IO-9 ) *12+IBLOCK+4) 

CALL MULT (Tl, QUl, DUN) 

CALL ADD (1 . , RI , -1 . , Rl / DUM) 

CALL MULT (DUM, Tl, DUO) 

CALL MULT (DUO, QS1, DUM) 

CALL ADD ( - 1 . , DUN , - 1 . , DUM , DUO ) 

CALL MULT (TIT, DUO, DUP) 

CALL MULT (T1,QS1, DUM) 

— CALL MULT (TIT, DUM, DUN) 

CALL ADD ( -1 . , DUN, 0 . , DUN, DUN) 

CALL JUXTV (DUP, DUN, DUO) 

555 CALL JUXTV (QU,QS, DUP) 

DO 85 IBLOCK =1,6 
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DO 84 JBLOCK =1,12 

ID= (APPEND* 6+IBLOCK-l) *NA+JBLOCK+BEAMl*12-8 
A ( ID) =DUO ( ( IBLOCK-1 ) *12+ JBLOCK+4 ) 

DO 85 JBLOCK =1,12 

ID= (APPEND* 6+ IBLOCK- 1 ) *NA+JBLOCK+ IBEAM* 12 -8 
A ( ID) =DUP ( (IBLOCK-1) *12+ JBLOCK+4 ) 

APPEND=APPEND+ 1 
CONTINUE 
CONTINUE 
EM=MASS ( OBODY+ 4 ) 

CALL SET ( INRTO ,3,3) 

DO 24 IBLOCK =1,9 

INRTO ( IBLOCK+4 ) =INRT ( (OBODY*3-3 ) *3+IBLOCK+4) 

CALL SET(PF,3, 12) 

CALL SET (PM, 3, 12) 

DO 33 IBLOCK = 1,36 

PF (IBLOCK+4) =P ( (IBEAM*12-6) *12+IBLOCK+4 ) 

PM ( IBLOCK+ 4 ) = P ( (IBEAM*12-3 ) *12 + IBLOCK+4 ) 

W2=W*W 

W2IN=1./W2 

CALL MULT ( RO , TBEAM , DUM ) 

CALL MULT (TIT, DUM, TRT) 

FIRST BLOCK, FORCE EQUATIONS, OUTBOARD 
IF (EM-999999999 .) 96, 97, 96 
AD=ENREF 
AE=0 . 

GO TO 98 
AD=ENREF * EM 
AE=W2 IN 
CONTINUE 

CALL ADD (AD, QU, AE, PF , DUO) 

CALL MULT (TRT, QS, DUN) 

CALL ADD ( 1 . , DUO , +AD , DUN , DUO ) 

SECOND BLOCK, MOMENT EQUATIONS, OUTBOARD 
CALL MULT (RO, TBEAM, DUM) 

CALL MULT (DUM, PF, DUN) 

CALL MULT ( TBEAM , PM , DUM ) 

CALL ADD ( +W2 IN , DUM , +W2 IN , DUN , DUP ) 

CALL MULT (TBEAM, QS, DUM) 

CALL MULT (INRTO, DUM, DUN) 

CALL ADD ( 1 . , DUP , ENREF , DUN, DUP) 

CALL JUXTV ( DUO , DUP , DUO ) 

DO 7 IBLOCK =1,6 
DO 7 JBLOCK =1,12 

A ( (OBODY*6+IBLOCK-7 ) *NA+IBEAM*12+ JBLOCK-8 ) =DUO ( (IBLOCK-1) 
1*12+ JBLOCK+4) 

CONTINUE 

RETURN 

END 

SUBROUTINE BODFORM (CONFIG, NBEAM, BODREF , NBODY ) 

INTEGER CONFIG, BODREF 

DIMENSION CONFIG (5,3), BODREF (10,2) 

FORMAT (7110) 

JMAX=1 

DO 50 IBEAM = 1 , NBEAM 
J=CONFIG ( IBEAM, 2 ) 

J2=J*J 

IF (JMAX* JMAX- J2 ) 1 , 2 , 2 

JMAX2=J2 

CONTINUE 

IF ( J) 51,51, 52 

J=-J 

BODREF ( J, 1 ) = IBEAM 
BODREF (J, 2) =0 
J=CONFIG ( IBEAM, 3 ) 


52 



J2=J* J 

IF ( JMAX* JMAX- J2 ) 3 , 4 , 4 

3 JMAX2=J2 

4 CONTINUE 
IF ( J) 53,53, 54 

BODREF ( J , 1 ) = IBEAM 
BODREF ( J , 2 ) = 1 
54 CONTINUE 

50 CONTINUE 

AD=JMAX2 

NBODY=SQRT (AD) + . 01 
DO 63 I = 1 , NBODY 

63 PRINT 777,1, BODREF (1,1) , BODREF (1,2) 

RETURN 
END 

SUBROUTINE PFORM (W, L, Z , EIX, EIY, EAZ , EIP , MPL, IPL, FZ , AGKX, AGKY, 
1PF, PM) 

REAL L, MPL, IPL 
DIMENSION PF (40) , PM(40) 

707 FORMAT ( 7E14 . 5) 

W2=W*W 
W2IN=1. /W2 

IF (EIX-FZ*FZ* . 071) 8, 8, 9 

8 ARG=MPL/ (FZ-EIX*MPL*W2/AGKX) 

BXAB=W*SQRT (ARG) 

BXCD=18.42/L 

GO TO 10 

9 CONTINUE 

AD= . 5* (MPL*W2 /AGKX-FZ/EIX) 

AE=MPL*W2/EIX 
ARG=AD+ SQRT (AD*AD+AE) 

"" BXAB= SQRT ( ARG ) 

ARG=-AD+SQRT (AD*AD+AE) 

BXCD=SQRT (ARG) 

10 CONTINUE 

IF (EIY-FZ*FZ* . 071) 11,11,12 

11 ARG=MPL/ (FZ-EIY*MPL*W2/AGKY) 

BYAB=W* SQRT (ARG) 

BYCD=18.42/L 

GO TO 13 

12 CONTINUE 

AD= . 5* (MPL*W2 /AGKY-FZ/EIY) 

AE=MPL*W2/EIY 

ARG= AD+ SQRT (AD*AD+AE ) 

BYAB=SQRT (ARG) 

ARG= - AD+ SQRT (AD*AD+AE) 

BYCD= SQRT (ARG) 

13 CONTINUE 

BZ= SQRT (MPL /EAZ) *W 
BP=SQRT (IPL/EIP) *W 
C FORCE MATRIX FOR BEAM 

CALL SET ( PF , 3 , 12 ) 

BXAB2 =BXAB * BXAB 
BXCD2 =BXCD * BXCD 
BXAB3 =BXAB2 *BXAB 
BXCD3 =BXCD2 *BXCD 
BYAB2=BYAB*BYAB 
BYCD2=BYCD*BYCD 
BYAB3 =BYAB2 *BYAB 
— BYCD3 =BYCD2 *BYCD 

IF (Z) 2 , 3 , 2 

PF ( 5 ) =+BXAB3 *EIX+BXAB*FZ 
PF (6) =0 . 

PF (7 ) =-BXCD3 *EIX-BXCD*FZ 
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PF ( 8 ) = 0 . 

PF (21) =+BYAB3 *EIY+BYAB*FZ 
PF (22 ) =0 . 

PF (23) =-BYCD3*EIY-BYCD*FZ 
PF (24) =0 . 

PF ( 3 7 ) =BZ * EAZ 
PF (38) =0 . 

GO TO 4 

PF (5) =-BXAB3*EIX*COS (BXAB*L) -BXAB*FZ*COS (BXAB*L) 

PF ( 6 ) = BXAB3*EIX*SIN(BXAB*L) +BXAB*FZ*SIN (BXAB*L) 

PF (7 ) = BXCD3 *EIX*COSH (BXCD*L) -BXCD*FZ*COSH (BXCD*L) 

PF ( 8 ) = BXCD3 *EIX*SINH (BXCD*L) -BXCD*FZ*SINH (BXCD*L) 

PF (21) = -BYAB3 *EIY*COS (BYAB*L) -BYAB*FZ*COS (BYAB*L) 

PF (22) = BYAB3 *EIY*SIN (BYAB*L) +BYAB*FZ*SIN (BYAB*L) 

PF (23 ) = BYCD3 *EIY*COSH (BYCD*L) -BYCD*FZ*COSH (BYCD*L) 

PF (24)= BYCD3 *EIY*SINH (BYCD*L) -BYCD*FZ*SINH (BYCD*L) 

PF (37) =-BZ*EAZ*COS (BZ*L) 

PF (38) =BZ*EAZ*SIN (BZ*L) 

4 CONTINUE 

C MOMENT MATRIX FOR BEAM 

CALL SET (PM, 3,12) 

IF(Z) 5, 6, 5 

6 PM (17) =0 . 

PM ( 18 ) =+BXAB2 *EIX+FZ 
PM ( 1 9 ) = 0 . 

PM(20) =-BXCD2*EIX+FZ 
PM(9) =0 . 

PM ( 10 ) =-BYAB2 *EIY-FZ 
PM (11) =0 . 

PM ( 12 ) =+BYCD2 *EIY-FZ 
PM (39) =BP*EIP 
PM (40) =0 . 

GO TO 7 , 

"5- PM(17) =-BXAB2*EIX*SIN (BXAB*L) -FZ*SIN (BXAB*L) 

PM (18) =-BXAB2 *EIX*COS (BXAB*L) -FZ*COS (BXAB*L) 

PM (19) =BXCD2 *EIX*SINH (BXCD*L) -FZ*SINH (BXCD*L) 

PM(20) =BXCD2*EIX*COSH(BXCD*L) -FZ*COSH (BXCD*L) 

PM(9) =BYAB2*EIY*SIN(BYAB*L) +FZ*SIN (BYAB*L) 

PM (10) =BYAB2*EIY*COS (BYAB*L) +FZ*COS (BYAB*L) 

PM(ll) =-BYCD2*EIY*SINH (BYCD*L) +FZ*SINH (BYCD*L) 

PM(12 ) =-BYCD2*EIY*COSH (BYCD*L) +FZ*COSH (BYCD*L) 

PM(39) =-BP*EIP*COS (BP*L) 

PM (40) =BP*EIP*SIN (BP*L) 

7 CONTINUE 
RETURN 
END 

SUBROUTINE QFORM (W,L,Z,EIX,EIY, EAZ , EIP , MPL , I PL , FZ , AGKX , AGKY , QU , QS ) 
REAL L , MPL , IPL 
DIMENSION QU(40) ,QS (40) 

W2=W*W 

IF (EIX-FZ*FZ* .071)1,1,2 

1 ARG=MPL/ (FZ-EIX*MPL*W2/AGKX) 

BXAB=W* SQRT (ARG) 

BXCD=18 .42/L 
GO. TO 3 

2 CONTINUE 

AD= - 5* (MPL*W2 /AGKX-FZ/EIX) 

AE=MPL*W2/EIX 
ARG=AD+ SQRT (AD*AD+AE) 

BXAB= SQRT (ARG) 

ARG=-AD+SQRT (AD*AD+AE) 

BXCD=SQRT (ARG) 

3 CONTINUE 

IF (EIY-FZ*FZ* .071)4,4,5 

4 ARG=MPL/ (FZ-EIY*MPL*W2 /AGKY) 



BYAB=W* SQRT ( ARG ) 

BYCD=18 . 42 /L 
GO TO 6 
CONTINUE 

AD= . 5* (MPL*W2 / AGKY-FZ/EIY) 
AE=MPL*W2/EIY 
ARG=AD+SQRT (AD*AD+AE) 
BYAB=SQRT (ARG) 
ARG=-AD+SQRT (AD*AD+AE) 
BYCD=SQRT (ARG) 

CONTINUE 

BZ=SQRT (MPL/EAZ ) *W 
BP=SQRT ( IPL/EIP) *W 
AD=1 . 

LINEAR DEFLECTION MATRIX 
CALL SET (QU, 3,12) 

IF(Z) 8, 9, 8 
QU (5) =0 . 

QU ( 6 ) =1 . 

QU ( 7 ) = 0 . 

QU (8) =1 . 

QU (21) =0 . 

QU (22 ) =1 . 

QU ( 2 3 ) = 0 . 

QU ( 2 4 ) = 1 . 

QU ( 3 7 ) = 0 . 

QU (38) =1 . 

GO TO 10 
CONTINUE 

QU (5) =SIN (BXAB*L) 

QU ( 6 ) =COS ( BXAB * L ) 
QU(7)=SINH(BXCD*L) 

QU ( 8 ) =COSH (BXCD*L) 

QU (21) =SIN (BYAB*L) 

QU (22 ) =COS (BYAB*L) 

QU{23 ) =SINH (BYCD*L) 

QU ( 2 4 ) =COSH ( BYCD*L) 

QU (37 ) =SIN (BZ*L) 

QU (38) =COS (BZ*L) 

CALL ADD(AD,QU, 0. , QU,QU) 
ANGULAR DEFLECTION MATRIX 
CALL SET (QS, 3, 12) 

IF (Z) 11 , 12 , 11 
QS ( 17 ) =BXAB 
QS (18) =0 . 

QS (19) =BXCD 
QS (20 ) =0 . 

QS ( 9 ) =-BYAB 
QS (10) =0 . 

QS (11) =-BYCD 
QS ( 12 ) =0 . 

QS (39 ) =0 . 

QS (40 ) =-l . 

GO TO 13 
CONTINUE 

QS (17) =BXAB*COS (BXAB*L) 
QS(18)=-BXAB*SIN(BXAB*L) 
QS (19) =BXCD*COSH(BXCD*L) 
QS (20 ) =BXCD*SINH (BXCD*L) 
QS (9) =-BYAB*COS (BYAB*L) 

QS (10) =BYAB*SIN(BYAB*L) 

QS (11) =-BYCD*COSH (BYCD*L) 
QS (12) =-BYCD*SINH (BYCD*L) 
QS (39 ) =-SIN (BP*L) 

QS (40 ) =-COS (BP*L) 

CALL ADD ( - 1 . , QS , 0 . , QS , QS ) 



RETURN 

END 


r 



SUBROUTINE PQFORM(W, L, E, P, Q, DUM, DUN) 

REAL L, LI , MPL, IPL 

DIMENSION L(9),E(5),P(5),Q(5), DUM ( 5 ) , DUN ( 5 ) 

NBEAM=L ( 1 ) + . 0 1 

DO 1 IBEAM = 1 , NBEAM 

LI=L ( IBEAM+4 ) 

EIX=E ( ( IBEAM- 1) *9+5) 

EIY=E( (IBEAM- 1) *9+6) 

EAZ=E ( ( IBEAM- 1) *9+7) 

EIP=E( (IBEAM-1) *9+8) 

MPL=E ( (IBEAM-1) *9+9) 

IPL=E( (IBEAM-1) *9+10) 

FZ=E ( (IBEAM-1) *9+11) 

AGKX=E( (IBEAM-1) *9+12) 

CALL PFORM^W^LI ! 0 . , EIX , EIY , EAZ , EIP , MPL , IPL , FZ , AGKX , AGKY , DUM 
1 , DUN) 

IF ( IBEAM- 1 ) 2 , 3 , 2 
CALL MAKE (P, DUM) 

GO TO 4 

CALL JUXTV ( P , DUM , P ) 

CALL JUXTV ( P , DUN , P ) 

CALL QFORM (W, LI , 0 . , EIX, EIY, EAZ , EIP , MPL , IPL, FZ , AGKX, 

1 AGKY , DUM , DUN ) 

IF (IBEAM-1) 5, 6, 5 
CALL MAKE (Q, DUM) 

GO TO 7 

CALL JUXTV ( Q , DUM , Q ) 

CALL JUXTV (Q, DUN, Q) 

CALL PFORM (W, LI , LI , EIX, EIY , EAZ , EIP , MPL, IPL, FZ, AGKX, 

1 AGKY, DUM, DUN) 

CALL JUXTV ( P , DUM , P ) 

CALL JUXTV ( P , DUN , P ) 

CALL QFORM (W, LI, LI, EIX, EIY, EAZ , EIP, MPL, IPL, FZ , AGKX, 

1 AGKY , DUM , DUN ) 

CALL JUXTV (Q, DUM, Q) 

CALL JUXTV ( Q , DUN , Q ) 

RETURN 

END 


SUBROUTINE WSEARCH ( W , DW , A , DETNEW , DETOLD , WI ) 
DIMENSION A ( 5 ) , DETOLD ( 5 ) , DETNEW ( 5 ) , WI ( 5 ) 
FORMAT (7E12 . 5) 

NA=A ( 1 ) 

A7 7 = A (77) 

CALL MAKE (DETOLD, DETNEW) 

DETO=l . 

DO 3 I =1 , NA 
DETO=DETO* DETOLD ( 1+4 ) 

CALL UPPER (A, DETNEW) 

DETN=1 . 

DO 1 1=1 ,NA 
DETN=DETN* DETNEW ( 1+4 ) 

FREQ=W/6 .2831853 

PRINT 707, FREQ , A7 7 , DETN 

IF (DETO*DETN) 4,2,2 

WROOT= W - DW- DETO * DW / (DETN-DETO) 

NW=WI (1) +1.001 
WI ( 1 ) =NW 
WI (NW+4 ) =WROOT 
CALL SPIT (WI , 3H WI) 

RETURN 

END 



SUBROUTINE UPPER (A, DETA) 

C GERNERATES THE DETERMINANT, DETA, OF MATRIX A 

DIMENSION A ( 5 ) , DETA ( 5 ) 

'1 FORMAT (2 El 5 . 7 ) 

_ N = A ( 1 ) + .01 

NMl = N - 1 
NPl = N + 1 
DO 5 K = 1 , NMl 
KROW=K 

AKK=A (K*N+K-N+4 ) 

IF (AKK) 1,2,1 

2 KROW=KROW + 1 
IF(KROW-N) 7,7, 5 

7 AKK=A(KROW*N+K-N+4) 

IF (AKK) 3,2,3 

3 DO 4 J=K, N 
AD=A (K*N+ J-N+4 ) 

A (K*N+J-N+4 ) =A (KROW*N+ J-N+4 ) 

4 A(KROW*N+J-N+4) =AD 

1 AKKl = 1 . /AKK 
K1 = K + 1 

DO 6 KI = Kl , N 

AKIK = A (KI*N + K - N+4)*AKKl 
A(KI*N + K - N + 4 ) =0 . 

DO 6 L = K1,N 

KIL = KI*N + L - N + 4 

A(KIL) = A (KIL) - AKIK*A (K*N + L - N + 4) 

6 CONTINUE 

5 CONTINUE 
AD=1 . 

CALL SET (DETA, N, 1) 

DO 16 I = 1 , N 

— DETA (1+4) =A(I*N + I - N + 4) 

EYE=I 

AD=AD*DETA ( 1+4 ) 

IF (AD) 26,27,26 
27 PRINT 101, AD, EYE 
AD=1 . 

2 6 CONTINUE 

16 CONTINUE 

RETURN 

END 

c 

SUBROUTINE DIAG (A, SHAPE) 

DIMENSION A (5 ), SHAPE (5) 

123 FORMAT (110, 5E12 .4) 

NA=A ( 1 ) + . 0 1 

NAM1=NA-1 

DO 1 I = 1 , NAM1 

IP1=NA-I+1 

AD=0 . 

DO 2 J = IP1 , NA 

2 AD=AD-SHAPE ( J+4 ) *A ( (NA-I-1) *NA+J+4) 

SHAPE (NA-I+4 ) =AD/A ( (NA-I-1) *NA+NA-I+4) 

1 CONTINUE 

RETURN 
END 
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ABSTRACT 

During the past three decades the finite element method matured gradually and dominated 
almost exclusively all the engineering applications. But the practice of generating complex finite 
element dynamic models of aerospace structures has revealed a number of shortcomings. First, 
the high dimensionality of the models requires an order-reduction process before a control system 
can be designed. However, seemingly unimportant modes can be inadvertently eliminated which 
prove later to be significant to control system performance and stability. Second, use of finite 
element models for dynamic system identification is generally very difficult because of the 
extremely large amount of unknowns. Third, the structural damping is generally added ad hoc 
after generating an undamped model. Inaccurate damping results, especially for the modes which 
couple different types of motion. 

In contrast, distributed parameter models offer an alternative to finite element models for 
overall dynamic analysis and control synthesis for aerospace and aeronautical structures. First, 
the modal order does not have to be reduced prior to the inclusion of control system dynamics, 
which eliminates the risk involved with modal truncation. Second, distributed parameter models 
inherently involve fewer parameters, thereby enabling more accurate parameter estimation using 
experimental data. Third, it is possible to include the damping in the basic model, thereby 
increasing the accuracy of the structural damping. Recently, distributed parameter models have 
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been made for some large space structures. Distributed parameter models may provide an 
efficient complementary methodology to the finite element approach. 

The software package PDEMOD was initialized by the late Lawrence W. Taylor, Jr. at 
NASA Langley Research Center during the middle of the 1980’s. The initial interest in the 
package PDEMOD was to model the structural dynamics of general spacecraft configurations by 
using the distributed parameter approach, which consists of a three-dimensional network of 
flexible beams and rigid bodies. The building blocks from which three-dimensional configurations 
can be constructed consist of (1) beams, which have bending in two directions, torsion, and 
elongation degrees of freedom, and (2) rigid bodies, which are connected by any network of beam 
elements. The full six degrees of freedom are allowed at either end of the beam. Rigid bodies can 
be attached to the beam at any angle or body location. The modified Bemoulli-Euler beam 
equation is used to represent the bending, the wave equations for torsion and elongation. 

A system of partial differential equations (PDEs) is formulated and connected at the 
elements’ boundaries based on the compatibility conditions. The equations of motion for any 
number of rigid bodies are written in the frequency domain and in terms of the coefficients of the 
sinusoidal and hyperbolic functions which comprise the mode shapes. The force and moment 
vectors for both ends of a single beam element can be described in terms of the spatial derivatives 
of the solutions of the corresponding PDE’s. Distributed parameter models can therefore be 
generated for any three-dimensional configurations describable by PDEs joined at their 
boundaries. Because of Mr. Taylor’s sudden demise, it becomes an urgent task to summarize and 
sift his research achievements and make them available to the other researchers. We have 
continued his work and recovered the basic functions of this package which may provide an 
opportunity to more researchers to apply distributed parameter modeling techniques to a variety 
of aerospace structures. The verification of the code has been conducted by comparing the results 
with those examples for which the exact theoretical solutions can be obtained. 

By investigating the potential of the distributed parameter modeling technique, we expect 
that the PDEMOD may be further developed in the following aspects: (1) structural dynamics, 
modal frequencies and mode shapes; (2) parameter estimation of modal characteristics; (3) 
structural damping; (4) control system dynamics; and (5) design optimization. In its present 
stage, however, only the first of the functions has been completed and included in the package. 
Meanwhile, a massive effort is being conducted which is expected to modify the mathematical 
model and global system generating procedure to develop the methodology for control analysis 
purpose. Instead of using the coefficients of the solution functions, the transfer matrix may be 
used finally, which provides a much more convenient way to describe the state-vector transition 
from one point of the structure to the other. All these research outputs are planned to be contents 
of the modified version of the package. 


INTRODUCTION 

During the past three decades the finite element method had become extremely popular in 
a very wide field of engineering. As early as the 1960’s, the aircraft industry had developed in- 
house finite element programs. During the 1970’s, some general-purpose finite element programs 
such as NASTRAN were released for public use, bringing with them a significant technology base 
that led to development of numerous commercial finite element software systems. Later, the 
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various commercial packages were refined, and their technology base was expanded. NASA s 
research in computational structures technology (CST)[i] is helping to develop the finite element 
analysis to a new stage that is capable of using a general automated unstructured grid generation 
to discretize the aerodynamic field and the structure for thermal and structural analysis [2] . 

But such a treatment having a very large number of elements is usually very expensive and 
time consuming. In static finite element analyses, models having over 10,000 degrees of freedom 
(DOF) are not uncommon. In fact, a finite element model of a helicopter fuselage currently under 
design contains about 27,000 grid points and 149,000 DOF. It is economically unfeasible and 
normally unnecessary to conduct dynamic analysis with this many unknowns. Most of the 
dynamic models based on the finite element approach must resort to modal truncation techniques 
to reduce the number of unknowns prior to dynamic analysis. However, the control spillover 
resulting from the modal truncation may lead to degradation of control system performance, or 
even create instability [ 3 ]. 

The configuration of large space structures usually has the following common 
characteristics: extremely large dimension, light-weight design, high flexibility, rather uniform 
mass and stiffness distribution, low and closely spaced natural frequencies, and slight and/or 
improperly modeled damping. Structures with these characteristics are essentially distributed 
parameter systems by nature, and should be most accurately modeled as a continuously distributed 
mass and stiffness over the entire structural area rather than a sequence of finite mass elements 

coupled together as the finite element approach does. 

In contrast, distributed parameter modeling provides a very practical approach for overall 
dynamic analysis and control synthesis for aerospace and aeronautical structures. Recently, 
distributed parameter models have been made for some large space structures, such as, the 
Spacecraft Control Laboratory Experiment (SCOLE)[ 4 ], Solar Array Flight Experiment (SAFE)[ 5 j, 
Space Station, Freedom^, Low-power Atmospheric Compensation Experiment (LACE) satellite 
model [ 7 ], and the Aerospace Large Flexible Manipulatory. The fundamental advantage of using 
the distributed parameter approach is to decrease the number of unknowns significantly. In the 
preliminary design stage, the detailed structural design is often neglected; therefore, a simple but 
efficient global structural model may be more beneficial to weigh the trade-off between the 
performance and cost, between the structural penalty and control system consummation, etc. 
Distributed parameter models may provide an efficient complementary methodology to the finite 
element approach. 

The software package PDEMOD was initialized by the late Lawrence W. Taylor, Jr. at 
NASA Langley Research Center during the middle of the 1980’s. The first release of his work on 
PDEMOD package was in 1987 [9> i 0 ,ii]. The initial interest in the package PDEMOD was to model 
the structural dynamics of general spacecraft configurations by using the distributed parameter 
approach, which consists of three-dimensional network of flexible beams and rigid bodies. The 
building blocks from which three-dimensional configurations can be constructed consist of (1) 
beams, which have bending in two directions, torsion, and elongation degrees of freedom, and (2) 
rigid bodies, which are connected by any network of beam elements. The full six degrees of 
freedom are allowed at either end of the beam. Rigid bodies can be attached to the beam at any 
angle or body location. The modified Bemoulli-Euler beam equation is used to represent the 
bending, the wave equations for torsion and elongation. 

A system of partial differential equations (PDEs) is formulated and connected at the 
elements’ boundaries based on the compatibility conditions. The equations of motion for any 
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number of rigid bodies are written in the frequency domain and in terms of the coefficients of the 
sinusoidal and hyperbolic functions which comprise the mode shapes. The force and moment 
vectors for both ends of a single beam element can be described in terms of the spatial derivatives 
of the solutions of the corresponding PDE’s. Distributed parameter models can therefore be 
generated for any three-dimensional configurations describable by PDEs joined at their 
boundaries. The manual labor of generating such models is therefore avoided. 

Because of Mr. Taylor’s sudden demise, it becomes an urgent task to summarize and sift 
his research achievements and make them available to the other researchers. We continued his 
work and recovered the basic functions of this package which may provide an opportunity to 
more researchers to apply distributed parameter modeling techniques to a variety of aerospace 
structures. The verification of the code has been conducted by comparing the results with those 
examples for which the exact theoretical solutions can be obtained, for instance, a simple beam 
with various boundary conditions, etc. 

By investigating the potential of the distributed parameter modeling technique, we expect 
that the PDEMOD may be further developed in the following aspects: (1) structural dynamics, 
modal frequencies and mode shapes; (2) parameter estimation of modal characteristics; (3) 
structural damping; (4) control system dynamics; and (5) design optimization. In its present 
stage, however, only the first of the functions has been completed and included in the package. 
Meanwhile, a massive effort is being conducted which is expected to modify the mathematical 
model and global system generating procedure]^] to develop the methodology for control 
analysis purpose] 13 , 14 ]. Instead of using the coefficients of the solution functions, the transfer 
matrix may be used finally, which provides a much more convenient way to describe the state- 
vector transition from one point of the structure to the other. All these research outputs are 
planned to be contents of the modified version of the package. 


PARTIAL DIFFERENTIAL EQUATIONS (PDEs) 

A complex large space structure can be decomposed into simple pieces. A network of 
distributed parameter elements and the attached rigid bodies are connected to represent the 
structural dynamics of the complex flexible spacecraft. Each flexible beam element exhibits lateral 
bending u and v in two axes, axial deformation w, and torsion \|/, as shown in Fig. 1 , which can be 
independently described by a variety of PDEsfisj. 



Figure 1 A Beam Element 


The bending behavior of a beam element can be described by a modified Bemoulli-Euler 
beam equation which includes Euler bending stiffness, Timoshenko shear, and axial-force stiffness 
for lateral deflections in the x-z and y-z planes. The PDE for bending in the x-z plane is 
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mu + El x u" + GAu' + F 0 u* = q x (z,t) 


( 1 ) 


The corresponding equation for bending in the y-z plane is 

mv + EI y v"+GAv'+F 0 v' = q y (z,t ) (2) 

The axial and torsional dynamics are represented by wave equations. 


mw- EAw’ = F x (z,t) (3) 

and 

J Y \j/ - Glyiif' = M t (z, t) (4) 

respectively. The PDEs provide the relationships between the modal frequency and the 
eigenvalues for the mode shape functions. The Euler and wave equations can be solved for the 
zero damping cases to produce the following relationships between the modal frequency oj and 
the eigenvalue P[i6], For bending in the x-z plane, 


B 2 =±- 
2 


mco 




GA El 


* J 


1 

*•' 1 ' 


r . 


mco 


mco 


GA El x I EI X 


For bending in the y-z plane, 


1 

B =±- 
2 


f 


mco 


GA El 


y ) 



mco 


V 


GA El 


y ) 


mco 

~EL 


For elongation and torsion, we have 
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( 6 ) 
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( 8 ) 


MODE SHAPE FUNCTIONS 

The solutions of these partial differential equations for zero damping produce the 
sinusoidal and hyperbolic spatial equations which comprise the mode shape functions. For the 
case that F 0 =0, the bending mode shape in the x-z plane is. 
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( 9 ) 


u — A x sin p x z + B x cos fi x z + C x sinh f5 x z + D x cosh /} x z 
Similarly, for the bending in the x-z plane the mode shape function is, 

v = A y sin fi y z + B y cos (5 y z + C y sinh j5 y z + D y cosh /3 y z ( 1 0) 

The undamped mode shape function for elongation along the z-axis is 

w = A z sin (3 z z + B x cosP z z (11) 

The undamped mode shape function for torsion about the z-axis is 

y= A y sin J3 y z +B y cos f5 y z (12) 

These undamped mode shape functions are expected to be good approximations to the 
exact solutions for low level of damping. The mode shape of the entire configuration consists of 
all these functions, repeated for each beam element. Because of bending in two directions, 
elongation, and torsion, a total of 12 coefficients are needed for each beam element. A vector of 
the coefficients of these sinusoidal and hyperbolic functions is defined as the mode shape 
parameter vector, 

{0} = [A,B,C,D'A t B,C t D,A,B,A r B r } T (13) 


The translational deflection vector is defined as. 
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(14) 


the non-zero elements of the matrix [Q u (z)] are as follows, 

Q " = sin p x z 0' u 2 = cos (l x z Q™ = sinh&z 

Ql s = sin P y z Ql 6 = cos (5 y z Q? = sinh fi y z 

Ql* = sin j5 z z = cos 

The angular deflection vector is defined as, 


Ql A =coshj3 x z 
Q 2t = cosh p y z 
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(15) 


v' =[e,(*)]{0) 

V. 

where, 

Q\ x Q' 2 Q] 3 Q' 4 0 0 0 0 0 0 0 0 

[&(z)] = 0 0 0 0 Q: 5 Q) 6 Q) 1 Q? 0 0 0 0 

|_ 0 0 0 0 0 0 0 000 Q 3U Ql' 12 . 

the non-zero elements of the matrix [Q s (z)] are as follows, 

Ql' = /3 x cosj3 x z Ql 2 =~P x sinp x z Q' s 3 = P x coshP x z Q\ 4 = /3 X sinh/3 x z 

Ql 5 = p y cos/^z QT=~P y sin (3 y z Q; 1 = p y cosh p y z 0” = /3, sinh/3 y z 

Ql'"= sin fa Ql u = cos P ¥ z 


FORCES AND MOMENTS 

Next, it is necessary to express the forces and moments at either end of the beam element 
in terms of the mode shape parameter vector {0}. The force vector is as 

>1 \EI x u~ 

{F} = ' F y ' = " EI y v m ' = [/V(r)]{0} (16) 

Fj EAw\ 

where, 

>» P' F 2 P F P F 0 0 0 0 0 0 0 0 

[ P F (z)]= 0 0 0 0 P 25 P 26 P F P F 0 0 0 0 

0 0 0 0 0 0 0 0 P 39 P 3 ' 10 0 0. 

the non-zero elements of the matrix [Pf(z)] are as follows, 

P? = ~PIEI X cos p x z P' F 2 = P\EI X sin P x z 

P' F 3 = p 3 x EI x cosh p x z P' F 4 = P 3 X EI X sinh/5 x z 

P? = ~Pl EI y C0S Py Z P F 6 = Pl EI y sin ft.* 

P; 7 = P 3 yEIy COSh P y Z Pf = P'yEIy S^ll PyZ 

P 39 = P : EA cos P t z P F ' 0 = ~P Z EA sin /3,z 

The moment vector can be expressed as 
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m) Ely 

{M} = j M y |= Ely = [p M (z)]{e] (17) 

m,\ [Giy, 

where, 

^ PS PS P'J 0 0 0 0 0 0 0 o' 

[/> M (z)]= 0 0 0 0 PS PS PS PS 0 0 0 0 

. 0 0 0 0 0 0 0 0 0 0 py py_ 

the non-zero elements of the matrix [Pm(z)] are as follows, 

Pm = s in fe PL 2 = -PIE*, cosj3 x z 

Pm = Pl EI * sinh P x z P" = p 2 x EI x cosh p x z 

Pm = ~K EI y sin Py z P m = ~K EI y cos / 3 y z 

Pm = sinh^r />” = p]EI y cosh P y z 

Pm" = P V GI V cos pyZ Py = -P y GI y sin PyZ 

It is also necessary to account for changes in axes from each beam to the body to which it 
is attached, and for points of attachment at some distance from the center of gravity of the body. 
The forces and moments that the beam i applies to the bodyj expressed in the body j’s coordinate 
system are, 

{y},=m,<F},=m,[y,] ( W}, (w) 

M,=[rJ,(tM), + W,{F} f )=[71 s ([p4 + [/(I ; ,[^] 1 ){e}, (19) 

where, the coordinate-transformation matrix from the beam i’s coordinate system to the body j’s 
coordinate system is 

cos(^,x,) cos(^ ; .,^.) cos(X ; ,z) 

lT\ji = cos(7 ; ,x,) cos (Y Jt y t ) cos(7 y ,z,) (20) 

_ cos(Z y , x,. ) cos(Z y , y - ) cos(Z y , z t ) . 

the eccentric matrix at the attachment point between the beam i and bodyj is 

0 —r z r r 

r z 0 -r x (21) 



8 



RIGID BODY MOTIONS 


A Newtonian or inertial frame of reference XoYoZo is used for the motions of all beam 
elements and rigid bodies. 



Figure 2 Rigid Body Motion 


Consider body j connecting several beams, of which beam k is taken into account as a 
datum. From Fig.2, we have 

at time t=to : R ao = R Co + r 0 

at time t=t: R a< = R aa +u,= R Ca +r 0 + u , 

On the other hand, R = R^ +r f , so, 

R c , = K, ~ 7 t =K +r 0 +u,-r t = R Co + u, +r 0 xu' t (22) 

where, u t and u' are the translational and angular deflection vectors at the attachment point 
between the body j and the beam k, respectively, and the vector difference r Q - r t has been 
expressed as the vector cross product of r a and i.e. Ar =r 0 —r t = r 0 x u', which can also be 
written in matrix form as 



Therefore, Eq.22 can be written in matrix form as 
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( 24 ) 


Differentiating Eq.24, we get the acceleration of the bodyj’s center of gravity (C.G.)> 

(25) 

The angular acceleration of the body j is simply expressed as 

{£},=m, *{«'},* ( 26 ) 


STRUCTURAL DYNAMIC EQUATIONS 

The equations of motion for the connected bodies and elements consist of blocks of terms, 
assembled in an order dictated by the body and beam indices. The mass times the acceleration of 
each body is related to the sum of forces caused by each beam element and each applied force. 

[m] j {i? CG } = ^{F} fi +£{/ } +{g} ( 2? ) 

J i m 


where, ^{F } is the sum of the i-beam forces acting on the body j; X{/} ;m 1S sum of the 

. } ' m . . . 
m-applied forces acting on the body j; {g} is the gravitational vector. It is similar for the moment 

equation. 




where, is the sum of the i-beam moments acting on the body j; 

is the sum of the moments acting on the body j caused by the i-beam forces {F} Jt ; ^{M} jn is 

n 

the sum of the n-applied moments acting on the bodyj; and ^ [ T] m; [ F] >m [ T] Jm {/ } jm is the sum of 


the moments caused by the m-applied forces acting on the bodyj. 

For the case of without applied forces and moments and neglecting gravity force, referring 
to Eqs.25 and 26, we derive the following two equations. From the force equation. 


From the moment equation. 


(29) 


(30) 
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Let us express both sides of Eqs.29 and 30 in terms of the mode shape parameter vector, then, 


-ffl’HflT'U&l +«Um„[&U<e>» = (3D 


Eqs.3 1 and 32 are the structural dynamic equations for the most general configurations. To 
demonstrate the overall procedure more clearly, let’s assume that there is only one beam element 
attached to the body j, that is, i=k=T. In this case, Eqs.3 1 and 32 will be simplified as 


[A F ]{d},={0} 

(33-1) 

3x12 12x1 


K]«}, ={oi 

3x12 12x1 

(33-1) 


where, 

M =m,[&],+ mja], +^[m];'[7-UiU 

Two equations in Eq.33 may be combined as 

U) {©}, = {0} (34) 

6x12 12x1 

where the system matrix [A] consists of the two block matrices [Ap] and [Am]. We can see that 
the number of equations is less than the number of unknown parameters. The difference is six. 
This is because there are six rigid-body degrees of freedom (d.o.f.s) for this particular one-body- 
one-beam system. We should have, therefore, six constraint equations to fix the six rigid-body 
d.o.f.s. For most common cases, six boundary conditions provide six constraint equations, which 
can also be expressed in terms of the mode shape parameter vector {0};, but the format of the 
equations depends on what the specific boundary conditions are. Superimposing the constraint 
equations into Eq.34, we obtain a new system dynamic equation 

={ 0 } (35) 

12x12 12x1 

which has a full-rank system matrix [Z] for any values of circular frequencies except the natural 
frequencies. Based on the condition that Det[A\ = 0, the natural frequencies of the system can be 
determined. 

For general configurations, the structural dynamic equation may become very 
complicated, but the procedure to generate the equation is the same as stated. In general, for a 
structural system consisting of J bodies and I beams, the structural dynamic equation has the form 
of 
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( 36 ) 


[A] {0} f ={0} 

(6*y)x(12*7) (J2*/)xl 


The difference between the number of unknowns and the number of equations is exactly equal to 
the number of rigid-body d.o.f.s. To fix the rigid-body d.o.fs, it is sometimes necessaiy to 
consider the compatibility conditions besides the boundary conditions. The following two special 
cases must be paid more attention. 

(1) For the case that a rigid body may have more than one beam element attached as shown in 
Fig.3 where two beams are attached to the body j, additional constraint equations must be added 
to the system equation, Eq.36, which appears as for this particular example in the figure, 

[A]{e} = { 0} (37) 

12x24 24x1 

It is clear that 12 extra equations are needed to fit the two bodies’ rigid-body d.o.f.s. While the 
boundary conditions provides six equations, the other six equations must be found from the 
compatibility conditions. 


ft 


beam ii 




— 

body j 



1 


[beam h 


Figure 3 A Rigid Body Attached by Two Beam Elements 


To account for the continuity in translational deflection at the two attachment points, the 
constraint equation must be satisfied. 




(38) 


The constraint equation for ensuring continuity in the angular deflection at the two attachment 
points is 


u %.<2 


(39) 


Expressing Eqs.38 and 39 in terms of the mode shape parameter vectors of the two beams, we 
have 


and 


(m,,,[e„]„ -[*]„, m, .„[&]>),, =(m w [al„ »[&],>>« 


(40) 

(41) 


Combining Eqs.40 and 41, we find the six additional equations. 
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(2) For the case that a beam connects two rigid bodies at its two ends as shown in Fig.4 
where beam i connects two bodies. More discussion must be addressed for this case. Apparently, 
it seems that the system equation, Eq.36, had a full-rank system matrix, since Eq.36 appears as for 
this particular example in the figure, 

[A]{d} i = { 0} ( 42 > 

12x12 I2xl 



Figure 4 A Beam Element Connecting Two Rigid Bodies 

But, it is wrong. The problem is that body ji and body j 2 connect to the beam i at different 
attachment points. Looking back to Eqs.14 to 17, and 31, we find that the components of the 
system matrix [A], such as [Q u (z)], [Q,(z)], [P F (z)], [P M (z)], are functions of the beam’s 
longitudinal coordinate z. Based upon the connection between body ji and beam i, a set of 
equations, which is the same as Eq.34, can be found, 

[x„(r, =0)]{6} i = {0} («-!) 

6X12 12x > 

Similar equation exists between body j 2 and beam i, 

[a, 2 (z, = I,)] {0) ■ = {0} (43-2) 

6X12 12X1 

The two equations in Eq.43 are not, however, independent since the system matrices [Ap] and 
[Aj 2 ] are both related to the same beam element, beam i. It can be proven that these two matrices 
are related by a constant matrix [O] which can be derived from beam i’s PDEs, that is, 

[^ 2 (z, = A)]=[®lk,te=0)] < 44 > 

We can only, therefore, choose one set of equations from Eq.43, and the other six equations must 
be fixed by the corresponding boundary conditions. 


VERIFICATION EXAMPLES 

1. EXAMPLE 1: Bending of a Cantilevered Beam. 
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Length of the Beam L=130.0in. 

Bending Stiffness EI=4xl0 7 Lb. in 2 . 

Mass per Length m=0. 09556 Lb.sec 2 /in. 


1.2106a„. The coefficients a n 


[W 

Theoretical formula for circular natural frequency: 0)„ = = 

and the theoretical natural frequencies and the corresponding results from PDEMOD are listed in 
Table 1. 


Table 1 Results of Example 1 


No. of Modes 

a* 

Natural Frequency 

Theoretical Value 

PDEMOD 

©n 

fn, HZ. 

fn, HZ. 

1 i 

3.52 

4.2613 

0.6782 

0.6775 

2 

22.0 

26.6332 

4.2388 

4.246 

3 

61.7 

74.6940 

11.8879 

11.890 

4 

121.0 

146.4826 

23.3134 

22.570 

5 

200.0 

242.1200 

38.5346 

38.490 


2. EXAMPLE 2: Bending of a Clamped-Clamped Beam. 


Length of the Beam L=130.0 in. 

Bending Stiffness EI=4xl0 7 Lb. in 2 . 

Mass per Length m=0. 09556 Lb.sec 2 /in. 

Theoretical formula for circular natural frequency: 

and the theoretical natural frequencies and the corresponding results from PDEMOD are listed in 
Table 2. 


[eT 

f: C0 " =a ^ = 


1.2106a„. The coefficients an, 


Table 2 Results of Example 2 


No. of Modes 

an 

Natural Frequency 

Theoretical Value 

PDEMOD 

(On 

fn, HZ. 

fn, HZ. 

i 

22.0 

26.6332 

4.2388 

4.31 

2 

61.7 

74.6940 

11.8879 

11.88 

3 

121.0 

146.4826 

23.3134 

23.29 

4 

200.0 

242.1200 

38.5346 

38.50 

5 

298.2 

361.0009 

57.4551 

57.50 


3. EXAMPLE 3: Bending of a Cantilevered Beam with a Tip-Mass M. 

Length of the Beam L=3 .077 ft. 

Bending Stiffness EI=1 75.9644 Lb. ft 2 . 
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Mass per Length m— 0.037 Lb. sec /ft (slug). 

Mass at the Tip M=0. 1538 Lb. sec 2 /ft (slug). 

Equivalent Stiffness of the Beam k=3EI/L 3 =l 8. 1202 Lb/ft. 

1 ^ 

Theoretical formula for the first circular natural frequency: ©, JJ+023m 


The theoretical 


natural frequency and the corresponding result from PDEMOD are listed in Table 3. 


Table 3 Results of Example 3 


No. of Modes 

Natural Frequency 

Theoretical Value 

PDEMOD 

co„ 

f n , Hz. 

fn,Hz. 

i \ 

10.566 

1.682 

1.736 


4. EXAMPLE 4: Torsion of a Cantilevered Beam. 

Length of the Beam L=130.0 in. 

Torsional Stiffiiess GI =4xl0 7 Lb.in 2 . 

T 2 

Polar Moment of Inertia per Length J,/L=0.2907 Lb. sec . 

Theoretical formula for the circular natural frequency: co„ = rm 
theoretical natural frequency and the corresponding result from PDEMOD are listed in Table 4. 


GI y 

f (J„ / L)L 2 


= 283.4745/?. The 


Table 4 Results of Example 4 


No. of Modes 

Natural Frequency 

Theoretical Value 

PDEMOD 

(On 

fn, HZ. 

fn,Hz. 

i 

283.4745 

45.1164 

45.10 

2 

566.9490 

90.2328 

90.19 


CONCLUDING REMARKS 

The computer software package PDEMOD is being developed to model the structural 
dynamics of general spacecraft configurations by using the distributed parameter approach. This 
paper provides the detailed description of the theoretical background used in the PDEMOD 
formulation. A complex large space structure is considered as an assembly of flexible beam 
elements and rigid bodies. Each flexible beam element is represented by four independent partial 
differential equations which exhibit lateral bending in two axes, axial deformation, and torsion. A 
system of partial differential equations is then formulated and connected at the elements 
boundaries based on the compatibility conditions. The equations of motion for any number of 
rigid bodies are written in the frequency domain and in terms of the mode shape parameter 
coefficients. The deflections, forces and moments for both ends of a single beam element can be 
described in terms of the spatial derivatives of the solutions of the corresponding PDE s, further 
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expressed in terms of the same set of mode shape parameter coefficients. Distributed parameter 
models can therefore be generated for any three-dimensional configurations describable by PDEs 
joined at their boundaries. The verification of the code has been conducted by comparing the 
results with those examples for which the exact theoretical solutions can be obtained. 

By investigating the potential of the distributed parameter modeling technique, we expect 
that the PDEMOD may be further developed in the following areas: (1) structural dynamics, 
modal frequencies and mode shapes; (2) parameter estimation of modal characteristics, (3) 
structural damping; (4) control system dynamics; and (5) design optimization. In present stage, 
however, only the first of these functions has been completed and included in the package. 
Meanwhile, a massive effort is being conducted which is expected to modify the mathematical 
model and global system generating procedure to develop methodology for the control analysis 
purpose. Instead of using the coefficients of the solution functions, the transfer matrix may finally 
be used, which provides a much more convenient way to describe the state-vector transition from 
one point of the structure to the other. All these research outputs are planned to be contents of 
the modified version of the package. It is also necessary to conduct additional testing to establish 
the accuracy of modeling different types of real space structures so as to move the approach from 
academic curiosity to a practical alternative for engineering design. 
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ABSTRACT 

The modeling, analysis and control of a small manipulator system are simplified by 
considering the links as rigid bodies. For large flexible manipulator systems, however, the 
flexibility of the links and the joint compliance must be considered. To describe the kinematics 
and the dynamic behavior of a flexible manipulator system, the common approach is to use 
Lagrange’s equations for both the rigid-body degrees of freedom (d.o.f.’s) and the dynamic 
deflection d.o.f.’s caused by the flexibility. The generalized coordinates are associated with the 
rigid-body d.o.f.’s of the links, and the modal coordinates associated with the flexibility d.o.f s. 
The consequence is that a set of highly-coupled and non-linear simultaneous partial differential 
equations is generated: These equations are so complex and lengthy that it is extremely difficult, 
if not impossible, to expand them manually even for a lower degree-of-freedom manipulator with 
a lower number of modes assumed. For the flexible manipulators with greater complexity, the 
dynamic analysis is literally forbidden by any practical manual symbolic derivations. The 
computer symbolic derivation of flexible manipulator dynamics was then suggested by several 
researchers. 

For simplifying the analytical process to a realistically acceptable extent, this paper 
conceives a new mathematical treatment for dynamic analysis of large flexible manipulator 
systems. The essence of the idea is to separate the kinematics and flexibility analyses as two 
independent but successive steps in a small time interval. Superposing the kinematic result and 
the flexibility effect, the summation is viewed as the initial conditions of the next instant motion, 
and the dynamic analysis succeeds to the next time interval. Repeating this process, the kinematic 
analysis accompanying flexibility effect is accomplished for a required time period. As usual, the 
kinematic analysis is based upon the rigid-body link assumption, and the Lagrange s equations are 
set up for these less amount of rigid-body d.o.f.’s, which are manually manipulative. The 
flexibility analysis for certain configuration of the manipulator system is conducted by using the 
distributed parameter system approach, along with the application of the transfer matrix method. 
Since an extremely complex analytical chore is resolved into two relatively simpler problems, the 
complexity of the dynamic analysis of large flexible manipulator systems is mathematically 
simplified. To demonstrate the applicability of the proposed methodology, an end-effector 
vibration suppression problem for a large manipulator system has been investigated. The 
manipulator system studied in this paper is a similitude of a NASA manipulator testbed for the 
research of the berthing operation of the Space Shuttle to the Space Station. The computational 
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results show that the proposed method is very effective for end-effector vibration suppression of a 
large flexible manipulator system. 


1. INTRODUCTION 

The modeling, analysis and control of a small manipulator system are simplified by 
considering the links as rigid bodies E i_4j. For large flexible manipulator systems, however, the 
flexibility of the links and the joint compliance must be considered. The limitations of the rigid 
link assumption in the formulation and analysis of large flexible manipulator dynamics were 
investigated extensively. Several formulations can be found in the robotics literature, such as, 
recursive or non-recursive Lagrangian assumed mode[5-7], generalized Newton-Euler method^], 
and Lagrangian using Raleigh-Ritz method^], etc. To describe the kinematics and the dynamic 
behavior of a flexible manipulator system, the common approach is to use Lagrange’s equations 
for both the rigid-body degrees of freedom (d.o.f.’s) and the dynamic deflection d.o.f.’s caused by 
the flexibility. The generalized coordinates are associated with the rigid-body d.o.f.’s of the links, 
and the modal coordinates associated with the flexibility d.o.f s. The consequence is that a set of 
highly-coupled and non-linear simultaneous partial differential equations is generated. These 
equations are so complex and lengthy that it is extremely difficult, if not impossible, to expand 
them manually even for a lower degree-of-ffeedom manipulator with a lower number of modes 
assumed. For the flexible manipulators with greater complexity, the dynamic analysis is literally 
forbidden by any. practical manual symbolic . derivations. The computer symbolic derivation of 
• ..flexible manipulator dynamics was later suggested by several researchers. Some of them wrote a 
■symbolic manipulation program^], some of them suggested using the MATHEMATICA 
commercial software package^]. The basic functions of the symbolic manipulation program may 
include symbolic simplification of polynomials and rational expressions, linearization of 
trigonometric functions, automated evaluation of the relative significance of terms and neglecting 
the less significant terms, and even symbolic integration and differentiation. The application of 
computer symbolic derivation techniques alleviates the difficulty in the dynamic analysis of large 
flexible manipulator systems. 

For simplifying the analytical process to a realistically acceptable extent, this paper 
conceives a new mathematical treatment for dynamic analysis of large flexible manipulator 
systems. The essence of the idea is to separate the kinematics and flexibility analyses as two 
independent but successive steps in a small time interval. Superposing the kinematic result and 
the flexibility effect, the summation is viewed as the initial conditions of the next instant motion, 
and the dynamic analysis succeeds to the next time interval. Repeating this process, the kinematic 
analysis accompanying flexibility effect is accomplished for a required time period. As usual, the 
kinematic analysis is based upon the rigid-body link assumption, and the Lagrange’s equations are 
set up for these less amount of rigid-body d.o.f.’s, which are manually manipulative. The 
flexibility analysis for certain configuration of the manipulator system is conducted by using the 
distributed parameter system approach, along with the application of the transfer matrix 
method[n]. Since an extremely complex analytical chore is resolved into two relatively simpler 
problems, the complexity of the dynamic analysis of large flexible manipulator systems is 
mathematically simplified. 

To demonstrate the applicability of the proposed methodology, an end-effector vibration 
suppression problem for a large manipulator system has been investigated. The manipulator 
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system studied in this paper is a similitude of a NASA manipulator testbed for the research of the 
berthing operation of the Space Shuttle to the Space Station^]. The system studied consists of 
two flexible links and three revolute joints which is assumed to be constrained in the vertical 
plane. There are only two rigid-body d.o.f.’s for this specific system. The two Lagrange’s 
equations are set up for the kinematic analysis, and the Runge-Kutta method is used to solve these 
non-linear partial differential equations numerically. The flexibility of the two links includes the 
lateral bending and axial elongation. The joint compliance is characterized by its torsional 
stiffness coefficient. The transfer matrices for the flexible arms and revolute joints have been 
constructed based on the partial differential equations. According to the compatibility conditions 
at the connecting points, the global system dynamic equation can be derived. From the 
corresponding boundary conditions, the characteristic equation for the global system is 
determined, from which the natural frequencies and mode shape functions can be found. 
Therefore, the transient response can be obtained. Joint moments are used as both displacement 
and control actuators. Control law computation proceeds in the frequency domain based on the 
pole-placement method^]. The computational results show that the proposed method is very 
effective for end-effector vibration suppression of a large flexible manipulator system. 

2. TWO-ARM FLEXIBLE MANIPULATOR SYSTEM 

The manipulator system (Figure 1) studied in this paper is a similitude of a NASA 
manipulator testbed for the research of the berthing operation of the space shuttle to the space 
station (Figure 2). This research testbed is planned to be the model of the berthing process 
constrained to move in the horizontal plane. Figure 2 illustrates the principal components of the 
facility. The Space Station Freedom (SSF) Mobility Base is an existing Marshall Space Flight 
Center (MSFC) Vehicle that has a mass of 2156.4kg. It represents a Space Station in the berthing 
operation. This vehicle is suspended on the MSFC flat floor facility using low flow-rate air 
bearings. The flexible appendage shown on the sketch will simulate solar panel disturbances. The 
vehicle has cold gas reaction jets to allow translational maneuvering. It also has a single gimbal 
for attitude control. The other vehicle, the Space Shuttle (SS) Mobility Base, is to be of similar 
construction and will be attached to the SSF Mobility Base with a flexible, two-link manipulator 
arm. The joints of the arms are driven by electric motors and are suspended by air bearings. 



Fig. 1 The Manipulator System Studied in the Paper 
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Fig. 2 NASA MSFC Manipulator Testbed 


The system studied consists of two flexible links and three revolute joints: shoulder joint 
A, elbow joint B and wrist joint C. In the paper, it is assumed that the base frame XoY 0 Z 0 is fixed 
on the Shuttle assumed as a rigid body, with Xo-axis along the joint A axis. The orientation of the 
Y 0 and Z 0 axes about the joint axis Xo is chosen such that the resultant base frame forms a right- 
hand coordinate system. The frame XiyiZi for link 1 is defined as follows: the xi-axis coincides 
with the joint axis of the joint A, the Zi-axis is along the axis of link 1. The frame x 2 y 2 Z 2 is fixed at 
the joint B with x 2 -axis coinciding with the joint B axis, rotating with link 2, and z 2 -axis being 
along the axis of the link 2. Since the dimension of the end-effector can not be in general 
comparable with the dimensions of the two links, the end-effector will be abstracted as a rigid 
body represented by a mass point as a whole at the joint C. 

Each joint of the manipulator arms is driven by an individual actuator. The control 
moments t A , Tb, and x c are also acting on the revolute joints A, B, and C, respectively. The joint 
compliance is characterized by its torsional stiffness coefficient. The corresponding input joint 
torques are transmitted through the arm linkage to the end-effector, where the resultant force and 
moment act upon the environment. The configuration of the corresponding rigid-body system of 
the one with flexibility can be specified by the two joint angles i f/i and I/O as shown in Figure 1 . 

3. KINEMATIC ANALYSIS UNDER RIGID-BODY LINK ASSUMPTION 


The common method for kinematic analysis under rigid-body link assumption in the 
robotics society is based on the solution of Lagrangian equation. 


dfdL' 

dt{dq ( j 



(3.1) 


where, Lagrangian function L is defined as L-T-U, T and LJ are the kinetic and potential energies 
of the system, respectively. Q t is the generalized force corresponding to the generalized 
coordinate qt. For the two-arm system discussed in this paper, two rigid-body d.o.f.’s are chosen 
as the two angles yr, and y/ 2 as shown in Figure 1, where points D and E are the mass centers of 
the two arms. The kinetic and potential energies for both arms can be expressed as follows, for 
arm 1: 


T x =\rny D +^I D y] 


U l = -m ] gL l sin\ir l 


(3.2) 
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For arm 2: 

T 2 =^m 2 v 2 B +^I E y/ 2 2 =^rn 2 [L\y/] + L x L 2 y/ x y/ 2 cos(y x - y/ 2 ) + ~L\y/\] 

U 2 = ]-m 2 g{lL x sin y/ x + L 2 sin y/ 2 ) 

Zr 

where, m h L t are mass and length of the ith arm (i=l or 2). The Lagrangian function is then, 


(3.3) 


l = t x + t 2 -u x -u 2 

1 . 1 


= -m x L] +~m 2 [L\ y] + 1, L 2 \j/ x \j/ 2 cos(^, - y 2 ) + - L\ y /\ ] 


+2 m 2 )L x sin y/ x +m 2 L 2 sin^J 


(3.4) 


Substituting Lagrangian function L into Eq.3.1 and taking corresponding derivatives, we derive 
the two Lagrangian equations: 


(-m. +m 2 )L\ i/>, +^-m 2 L x L 2 y/ 2 cos(^, - y 2 ) + -m 2 L x L 2 \j/\ sin(^, - y/ 2 ) 
3 2 z 

+K a y/ i ~K B y/ 2 +~ (m, + lm 2 )gL x cos <js i = r A - r B 


(3.5-1) 


— m 2 L x L 2 i// ] cos(y/ x - y/ 2 ) +~: rn 2 I ^¥ 2 ~ ~ sin(^, y / 2 ) 

2 s& 

1 

+K B y / 2 +~m 2 gL 2 cosy / 2 = r B -r c 

Isolating ”, and " 2 in each of the above two equations, Eq.3.5 can be reorganized as, 

^/n, + m 2 -\m 2 cos 2 (^, - + i w 2 ^i sin2(</, - y 2 ) + 2 m 2 L x L 2 y/ 2 sin(^y , - y/ 2 ) 

3 L. cos(^, - \f/ 2 ) 

+K A y , \-^+~ 


2Z,, 


]AT 5 ^ 2 + (ym, + m 2 )gL x cos^, - A m 2 gL x cosy / 2 cos (^, - (^ 2 ) 


= (L* -r s) - 


3i, (r 5 -r c )cos(^ 1 -y 2 ) 


2Z,, 


(3.5-2) 


(3-6-1) 


, 2( t wj, +m 2 ) 

^ co s(^-^)- 3cos( ^_^ 


L x L 2 ij r 2 + (y/”i + m 2 )Z^y'f tan((tq - ^ 2 ) + yW 2 Z, 1 L 2 ^ 2 sinf^, \y 2 ) 


+K aV\ ~ 


1 + 


2 (yffl[ + ^2)^1 
m 2^2 COS(^ - ^ 2 ) 


, ( 7 /W,+W2)^COS^ 2 

ft* + i<“ l +2m >>^ - co^.-r,) 


= (^ " r s)-‘ 


2(f m. +m 2 )I,(r s - r c ) 


m 2 Z. 2 cosO, - i/ 2 ) 


(3.6-2) 


5 



A set of two highly-coupled and non-linear simultaneous partial differential equations in Eq.3.6 
must be solved numerically. The well-known Runge-Kutta method is used. To do so, we define 


x, = y x it), x 2 = y/ x {t\ x 3 = y/ 2 it), x 4 = y/ 2 it) 

then, a set of four first-order simultaneous equations is generated, 

*1 = fi (') — *2 

x 2 = / 2 () = A x x\ sin2(x, -x 3 ) + £,x 4 sin(x, -x 3 ) + C,x l + £>,x 3 + £, cosx, 
+-F, cosx 3 cos(x, - x 3 ) + G, cos(x, - x 3 ) + H X 

X 3 — fz (') — ^4 

x 4 = / 4 (•) = A 2 x\ tan(x, - x 3 ) + B 2 x\ sin(x, - x 3 ) + C 2 x, + D 2 x 3 + E 2 cosx, 


+F, 


cosx. 


- + G. 


cos(x, -x 3 ) 2 cos(x, -X 3 ) 


+ H, 


where, 

3 A 

Eh. 

A, 


ilhkh.' C| ._£x. 

s " A, A, 


E x =- 


A, 

'jiPh + 2m 2 )gZ., ^ 4 m i8 L i 

A, 


E =■ 


G, =- 




M 


l ^A 2 

2A 2 A, 


H x =- 


*A~ X B 


1 3 2 2 

A, = [ , ffi, + m 2 ^ ^ cos (x, x 3 )]/., 


and 


(3.7) 


(3.8) 


A=~ 


1 2 
( T m, +m 2 )Z 1 


B, = - 


2A, 


C 5 =-— , d 2 =— [1 + 


2( 3 ot, +m 2 )L, 




(m, +2/?t ; )gL 1 


( 3 w, +w 2 )gZ 3 


2( 3 w,+w 2 )L 1 (r B -r c ) 

mjLjAj 


=~ ' F 2 = ‘ a * u 2“~ 

2A 2 A 2 

, 2( 3 /n,+m 2 ) 

A 2 =[ 2 m 2 cos(x, -x 3 )-- ; ~]L x L 2 

3cos(x,-x 3 ) 

The iteration formula of the fourth-order Runge-Kutta method is as follows^]: 


x i,j + 1 - x i.j + 2 +2k j3 + k i A ) 


i = 1 , 2, 3 and 4 


where 


(3.9) 


^i.l — fi ( X 1J > X 2J ’ X 3.j > X 4,j > 0 
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h , h h h h 

ki . 2 = y* (■'•1.J + ^1.1 > + ^ ^2,1 > *3.; + 2 ^3.1 > *4j + ^ ^ 4 -' ’ + ^ 

h , h h h h 

k it 3 = fi( x \,j + ~^k\ 2 > X 2.j + 2 kj.l^Z.j + 2 ^3.2 > X 4J + 2 4 2> ^ + 2^ 

k iA = /, (x, j + M, j ; x 2 ; + hk 2 3 ; x 3 ; . + M„ ; x 4J + M 4 3 ;/ + /») 


for the jth iteration, and h is the time interval. The accuracy of the method is in the order of h. 
The motion of the rigid-body manipulator system can now be solved by using the iteration 
formula, Eq.3.9. 


4. FLEXIBILITY ANALYSIS USING TRANSFER MATRIX METHOD 


In this paper, the inclusion of the flexibility of the manipulator arms is treated by the 
distributed parameter approach along with the transfer matrix method. The function of the 
transfer matrix is to relate the linear and angular deflections, forces and moments at one point in a 
structure to those at another point. The derivation of the transfer matrices for bending and 
elongation of a beam element is shown in Refs. 1 5 and 1 6, but is summarized here. The two 
flexible manipulator arms are represented by the Bemoulli-Euler equation and wave equation for 
their bending and elongation characteristics, respectively. The Bemoulli-Euler equation is in the 
form of 


tfUy 1 £U, 
dz 4 a 2 dt 2 


= 0 


(4.1) 


where, a y 2 =k y /m, and k y =EI is the bending stiffness and m is the mass per length of the beam. The 
elongation vibration is described by the wave equation 


d L u 1 1 i fu_ z 

~dF~~ai a 2 


(4.2) 


where, a z 2 =k z /m, and k z =EA is the axial stiffness. After separation of variables, Eqs.4.1 and 4.2 
can be expressed in the spacial domain as 

= 0 (43) 

0 (4.4) 

where, p y 4 = co 2 /a y 2 , p z = co/a z , and © is the circular natural frequency. The solutions to Eqs.4.3 and 
4.4 are 


(4.5) 


and 


U y (z) = A sin p y z + B cos p y z + C sinh p y z + D cosh fi y z 
U z (z ) = Msin/? 2 z + N cos fl 2 z 


(4.6) 
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where, A, B, C, D, M, N are the modal participation coefficients. If the state vector, {0}, is 
defined, then 

M = [ U,,U„r„F r ,F„M x ] T (4.7) 

where, U y and U z are the displacements along y- and z-axes, y/ x is the rotary angle of the beam 
sections about x-axis, F y and F z are the shear and tensile forces respectively, M x is the bending 
moment about x-axis. The state vectors at the two ends of a beam element are related by a matrix 
[<E>], that is, 


{<*£)) = [H{*(0)} 


(4.8) 


where, the transfer matrix [O] involving bending and elongation of a beam element as described in 
Ref. 16 is given in Eq.4.9, 


[o] = 


1 


A 

— (CO SfiyL 

+ cosh pyL) 

0 

1 

-Pyi-smpyL 
+ sinhy?yL) 

” ^ yfiy Py ^ 

+ sinh/3yL) 

0 


co sJ3 z L 


0 


2^ <Si " W 

+ sinhy9yL) 

0 

1 

“(co sB v L 
2 * 

-f cosh PyL) 


1 


v(- sin p v L 

2kyPy y 

+ s inhfiyL) 

0 

1 

Y (- cos p v L 

2 kyfip y 

+ cosh PyL) 

1 


(-COS PyL 


k z Pz 


'sin p z L 


~ kypy (- co spyL “ (cosp y L 


0 


0 


-k z P z sin p z L 


1 2 

kyPy ( COS P y ^ 


+ cosh/?^L) 
0 

1 


+ cosh pyL) 


^ k yP y ( sin PyL 
+ sinh/?^,L) 


+ cosh PyL) 
0 

1 

(sin PyL 


co sP z L 


UyPj, 

+ cosh PyL) 

0 

+ sinhPyL ) 

1 

~Py(r sin PyL 

+ sinh pyL) 

0 

1 


2 Py 

+ sinh PyL) 


^“(CO SPyL 
+ cosh PyL) 

(4.9) 

in which the elements at the 2nd and 4th rows and columns reflect the elongation, the others for 
the bending. Because the two arms are not in the same orientation, it is necessary to account for 
the alignment of the two arms. The relationship of the two local coordinate systems can be 
described by a coordinate transformation matrix [T 21 ], that is, 


1 y 2 
Iz., 


y 1 


-m; 


(4.10) 


where. 


fo,l = 


cosO^,) cosCy 2 ,z,) 
[_cos(z 2 ,.Vl ) cos(z 2 ,z,). 
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The transfer matrix of a revolute joint R has been derived in Ref. 11. A revolute joint R is 
abstracted as a massless torsional spring with spring constant k m through which two elements a 
and ‘b” are connected. The actuator fixed at the joint R will produce a control moment t r . The 
state vector {0} a at the end of element ‘k”is related to the state vector {0} b at the end of element 
“b” through the transfer matrix [O] of the joint R by 

(411) 

where, the transfer matrix [Or] of the joint R is given by 

"N 

N = 1 w 

1 


and the control-influence vector. 



0 0 


1 


V 2 + k Vt 


\T 


0 0-1 


(4.13) 


Applying the general expression, Eq.4.1 1, to the joints A and B, we find out that, 

{^.=o)} 0 =K]l^ z i =0 )}. + (^l^ 

and 

[6{z x = L,)}, = [O B ]{0(z 2 = 0)} 2 +{5 B }r B 

where, the subscripts 0, 1 and 2 stand for the Shuttle Base, arm 1 and arm 2, respectively. 

5. SYSTEM DYNAMIC EQUATIONS FOR A SPECIFIC CONFIGURATION 

For a specific configuration at an arbitrary time instant, the system dynamic equation was 
derived in detail in Ref. 1 1, that is 

fe=4)) 2 =[^]fe=0)} o+ [5]{r} (5-1) 

where, the system matrix [A]=[02][0 b] 1 [0 1 ][ < 3 ) a] the control-influence matrix ^ 

-[^[ObI^Bb}, -{Be}]; the control vector {x}=[ x A , tb, ^c] • 
The element transfer matrices [0.]’s and control influence vectors {B.}’s are associated with the 
two links and the three joints, respectively, recognized by the corresponding subscripts. 

The dynamic properties at a certain configuration without control actions are the inherent 
properties of the system, which are varied with the change in configuration when the manipulator 


(4.14) 

(4.15) 
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system is in motion. The method to derive these inherent dynamic properties is straightforward. 
Consider the boundary conditions (B.C.’s) of the system in Figure 1, at the fixed end (attachment 
point to the Shuttle Base), 

£/,„ (*, = 0) = U u (z, = 0) = ^ (2, = 0) = 0 (5.2) 


at the free end (end-effector), 

F yi (z 2 = L i ) = F Ii (z 2 = 4) = M Xi {z 2 = 4) = 0 


(5.3) 


Applying the B.C.’s, Eqs.5.2 and 5.3, to Eq.5.1 without control action, two equations can be 
derived, 


and 


= °) 

Wj ^,U=o) 
KU = o)J 


) = 


i U,( h = L,) 


(5.4) 


FyU, = 0 ) 

[4,]1 /;(*,= o) | =o 


(5.5) 


where, [A 12 ] and [A 22 ] are the block matrices of the matrix [A], The condition for Eq.5.5 having 
non-trivial solution is that 

DET[A 22 ]=0 (5.6) 

Eq.5.6 is the characteristic equation of the system from which the natural frequencies co's can be 
derived. After obtaining the natural frequencies from Eq.5.6, we can derive the mode shape 
functions from the equation below, 

[GJ{Q~0 (57) 

The detailed derivation of Eq.5.7 can be found in Ref. 16. The vector (Q consists of the modal 
participation coefficients for the beams 1 and 2 appearing in Eqs.4.5 and 4.6, that is,{Q=[Ai, Bi, 
Ci, Di, Mi, Ni, A 2 , B 2 , C 2 , D 2 , M 2 , N 2 ] t . Normalizing Eq.5.7 with N 2 =l, Eq.5.7 can be solved to 
obtain the modal participation coefficients for the two beams, thereby the mode shape functions 
can be obtained based on the solution equations, that is, Eqs.4.5 and 4.6, for the ith beam, 

U yi (z l )= A, sin P y z i +B i cos fi„z,+C, sinh fi^+D, cosh p y z t (5.8) 

and U h ( 2 , ) = si n z t + N t cos/? Z( z, (5 . 9) 

It is assumed that the control actions are related to the feedbacks of the nodal displacements and 
velocities, that is. 
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(5.10) 

(5.11) 

(5.12) 


r, = [*,]M*, = 0)}, = [*,, +**«,** + V-** + V.0-0,0]!^ = 0)}_ 

t b = [^b ]{^( z 2 = ^)} 2 = [^b, , k +k B} S,k B} +A: Bt 5,0,0,o] j#(z 2 — o) 

T C = [^cl{^( Z 2 = ^2 ) } 2 = [^"C, +M.^ + M.*Ci — "^2 ) } 2 

Inserting the control actions into Eq.5.1 and combining the similar terms, the closed-loop system 
dynamic equation follows. 


{0(z 2 =I 2 )} 2 =[^]^(z i =O)[ o 


(5.13) 


where, [Z] = [0 2 ][0 B ] ’[ojjo,] and 

[®,] = K] + { B , }[*,], [$„] = [o,] + {B, }[*,], [o,] = (i/] + {B c )[t c ])‘'[®,] . 

By applying the B.C.’s, Eqs.5.2 and 5.3, the closed-loop characteristic equation, DE1\A 22 \ = 0, 
can be derived, where, [4> 2 ] is a block matrix of the matrix \A\, from which the closed-loop poles 
can be found. 

6. SUPERPOSING RIGID-BODY KINEMATICS AND FLEXIBILITY EFFECT 

This, paper conceives a new mathematical treatment for dynamic analysis of large flexible 
manipulator systems. The essence of the idea is to separate the kinematics and flexibility analyses 
as two independent but successive steps in a small time interval. Section 3 of this paper gave the 
kinematic analysis assuming a rigid-body link based on the Lagrangian equation method, using the 
Runge-Kutta numerical approach. Sections 4 and 5 provided a method for system dynamic 
analysis due to flexibility at a specific configuration of the manipulator system. Compared with 
the macroscopic motion of the manipulator system, the motion resulted from flexibility is only a 
“microcosmic” motion. Only after a long-term effect is accumulated will the flexibility effect be 
significant. It allows us, therefore, to make the assumptions: at a certain configuration of the 
manipulator system, the deflections and the rates of deflection of the arms due to flexibility are 
small, and the elongation deformations are high-order infinitesimal so that they are neglected. 
Superposing the rigid-body motion and the motion due to flexibility, we have 

y } (z ]t t) = \j/^{t) + dy x {z x ,t) and y 2 (z 2 ,t) = V 2 (t) + dy 2 (z 2 ,t) (6.1) 

where, the small perturbation can be assumed as (cf. Fig.3), 

dy l (z l ,i) = —y } (z t ,t) and dy 2 (z 2 ,t) = ^-y 2 (z 2 ,t) (6.2) 

z, z 2 

Therefore, the two generalized coordinates defining the instantaneous motions of the two arms in 
the Lagrangian equation, Eq.3.5, would be 
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(6.3) 


Fig . 3 Small Perturbation of the Angles 
= t) and y/ 2 (t) = yt 2 (z 2 ,t)-— y 2 (z 2 ,t) 


Substituting Eq.6.3 into Eq.3.6 and neglecting high-order infinitesimals, we can express the 
Lagrangian equations in terms of and ~ 2 as follows, 

= A x y/ 2 x sin2(^, -y 2 )+B x y/\ sin(y7, -^ 2 ) + C,^, + D x y7 2 + E X cos^, 

+F X cosy/ 2 cos(^, - V 2 ) + G x cos (Jj/ X - y7 2 ) + H x + — y x (6.4-1) 


Yi = A 2 y/] tan (\f/ x - y/ 2 ) + B 2 if/\ sin(^, - y 2 ) + C 2 \j/ x +D 2 y/ 2 +E 2 cos y/ x 

cos u/ 7 „ 1 1 .. 

+F 2 — — — :zr— + G 2 — + H 2 + y 2 

cos(^, - y/ 2 ) cos(^, -y/ 2 ) z 2 


(6.4-2) 


Defining w x = y/ x , w 2 = y/ x , = y/ 2 , and w 4 = y/ 2 , a set of four first-order simultaneous 
equations is obtained which is similar to Eq.3.8, 


w x =w 2 


w 2 = A x w 2 sin2(H' 1 -w 3 ) + B x w 4 sin(w, -w 3 ) + C x w x + D x w 3 +E X cosw, 

1 .. 

+Ej cosw 3 cos(w, -w 3 ) + G, cos(w, -w 3 ) + H x + y x 


W 3 =w 4 


w 4 = A 2 m>\ tan(w, - w 3 ) + B 2 w] sin(w, - w 3 ) + C 2 w x + D 2 w 3 + E 2 cosw, 

+E, — 7~~ J z + Gj — r~ z + H 2 +—y 2 


cosw 3 [ _ 

cos(rv, -w 3 ) 2 cos(w, -w 3 ) 


(6.5) 


The only difference between Eq.6.5 and Eq.3.8 is the inclusion of y x and y 2 , which represents the 
effect of flexibility, and can be solved based on the formulation described in Sections 4 and 5 as 
long as the instantaneous configuration is specified and the motion at the end of previous time 
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interval is known. By using the Runge-Kutta procedure in Section 3, Eq.6.5 can be solved 
numerically. A solution is thus obtained which represents the superposition of rigid-body 
kinematics and flexibility effect. 


7. SIMULATION EXAMPLES 

As mentioned earlier, the manipulator system studied here is a similitude of a NASA 
MSFC manipulator testbed for the research of the berthing operation of the Space Shuttle to the 
Space Station. The mechanical properties of the system are listed in Table 1. 

Table 1 Mechanical Properties of the Two-arm Manipulator System 



Beam 1 

Beam 2 

Beam Length, L, (in.) 

60.0 

120.0 

Sectional Area, A, (in/) 

50.27 

12.57 

Second Moment of the Cross Section, I, (in. 4 ) 

107.0 

107.0 

Modulus of Elasticity, E, (psi) 

201.06 

12.57 

Mass per length, (slug/in.) 

0.1526 

0.0382 


The initial conditions of the system are assumedat y/ x =90" ,y/ 2 -0° ,y/ x - 0, and y/ 2 - 0 . 
First, the difference between the results of under rigid-body link assumption and with flexibility 
effects are noted by three examples as shown in Fig.4. These examples exhibit the motions of the 
manipulator system from the specified initial conditions under the actions of gravity force, inertia 
force, and the constant joint control moments Ta, Tb, and Xc as shown in each figure. The solid 
lines represent the motions under rigid-body link assumption, while the dashed lines represent the 
motions with flexibility effects. For clarity, the vibratory wave shapes for flexible beams were 
neglected in the figures. 

Next, several examples demonstrate the effectiveness for end-effector vibration 
suppression. The motion of the manipulator system is a continuous process, which stimulates 
vibration of the system at every time instant. The vibration suppression action continues 
throughout the whole process without interruption. The joint control moments will always 
change themselves based upon the control law given in Eqs.5.10 to 5.12. They play the roles of 
both actuating the manipulator system to fulfill a certain task and alleviating vibratory fluctuation 
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Fig.4 The Difference Between the Results Of Under Rigid-Body Link Assumption 

And With Flexibility Effects 

Figs. 5 to 7 demonstrate the effectiveness for the end-effector vibration suppression at 
instantaneous positions 1, 2, 3 as shown in Fig.4, with the initial joint control moments Ta= 1 ft- 
klb, Tb= 2 ft-klb, and tc=4 ft-klb. Both time histories without control (left) and with control 
(right) are shown in the figures for comparison. The vibratory motion of the end-effector is 
described in the second link’s coordinate system, X 2 y 2 Z 2 , as defined in Section 2. The upper two 
figures in Figs. 5 to 7 represent the results in y 2 -direction, the lower two in Z 2 -direction. The 
instantaneous vibration can be suppressed in about 0.3 second for all positions studied. 

8. CONCLUDING REMARKS 

A new mathematical treatment for dynamic analysis of large flexible manipulator systems 
is derived. An extremely complex analytical chore is resolved into two relatively simpler 
problems, the complexity of the dynamic analysis of large flexible manipulator systems is, 
therefore, mathematically simplified to a realistically acceptable extent for practical manual 
symbolic derivation of the equations of motion. Since the equation of motion is a set of highly- 
coupled and non-linear simultaneous partial differential equations, the Runge-Kutta numerical 
procedure has been used to solve the equations. As an example, the vibration suppression 
problem of a similitude of a NASA MSFC manipulator testbed has been investigated in the paper. 
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The computational results show that the proposed method is very effective for end-effector 
vibration suppression for a large flexible manipulator system. 
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Fig. 5 End-Effector Vibration Suppression at Instantaneous Position: y/ x = 103.2°, - 3.6° 
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Fig 6 End-Effector Vibration Suppression at Instantaneous Position: y / , = 121.8° ,y/ 2 = 10.2° 
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